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1  Introduction  and  Summary 


This  project  has  been  concerned  with  basic  research  on  analytical  methods  for  represen¬ 
tation  of  multiple  scattering  effects  in  heterogeneous  random  media.  We  have  examined 
several  methods  for  representation  and  computation  of  effective  parameter  models  for 
the  media  in  terms  of  statistical  characterizations  of  the  media.  These  methods  lead  to 
representation  for  the  effective  parameters,  e.g.,  the  effective  dielectric  constant  or  con¬ 
ductivity  in  a  random  medium.  Since  these  representations  may  be  difficult  to  evaluate 
in  specific  cases,  we  have  also  examined  current  research  on  the  derivation  of  bounds 
for  the  effective  parameters  as  a  means  of  quickly  providing  good  approximations  to 
the  behavior  of  the  medium.  Finally,  we  have  included  a  computational  technique 
for  the  evaluation  of  simplified  scattering  (and  absorption)  representations  in  terms 
of  transmission  and  reflection  coefficients  for  EM  radiation  encountering  a  region  of 
foliage. 

Our  primary  objective  has  been  analysis  of  the  interaction  of  EM  radiation  with 
microscopic  scattering  elements,  retaining  multiple  scattering  phenomena  in  the  final 
analytical  representation  of  the  physics  of  the  interaction.  Since  the  geometries  assumed 
by  fohage  as  a  scattering  media  are  random  and  may  include  several  different  types  of 
elementary  scatterers,  the  physics  is  extremely  complex.  Figure  1  shows  the  scattering 
cross  section  of  a  single  tree  in  the  frequency  range  of  interest  here.  It  is  apparent  from 
the  figure  that  the  infrastructure  of  the  tree  plays  a  significant  role  in  determining  the 
radar  cross  section. 

We  have  focussed  our  research  on  analytical  tools  capable  of  producing  systematic, 


computationally  feasible  approximations  to  the  multiple  scattering  process.  We  have 
concentrated  our  efforts  on  the  evaluation  of  effective  parameter  representations  of 
the  scattering  medium  which  retain  multiple  scattering  processes  in  the  definition  and 
computation  of  the  effective  parameter  models.  We  have  examined  three  different 
methods  for  this  purpose: 

1.  Methods  based  on  the  Transition  Matrix  formalism  which  represent  the  explicit 
interaction  of  elementary  scatterers  and  permit  the  systematic  construction  of 
approximations . 

2.  Bounds  for  the  effective  paramters  (dielectric  constant)  representing  the  medium 
based  on  the  generation  of  certain  “trial  fields”  and  the  use  of  geometry  dependent 
and  geometry  independent  formulations. 

3.  Homogenization  of  the  medium  via  an  asymptotic  analysis  of  Maxwell’s  equa¬ 
tions  in  differential  form  producing  an  approximation  to  the  total  fields  in  the 
heterogeneous  medium  in  terms  of  the  corresponding  fields  in  a  medium  with  a 
homogeneous  model  (e.g.,  dielectric  constant). 

Our  work  on  the  T-matrix  methods  is  summarized  in  Section  2.  It  is  important 
to  emphasize  that  the  T-matrix  methods,  though  formal  in  structure,  provide  simple 
models  for  multiple  scattering  phenomena  which  are  valid  for  a  wide  range  of  frequen¬ 
cies  (of  the  incident  radiation).  We  have  described  a  series  of  approximations  (ATA, 
CPA,  EMA,  SCA)  based  on  the  T-matrix  formalism.  We  have  also  shown  how  these 
approximations  can  be  developed  systematically  by  associating  the  T-matrix  formalism 


2 


with  a  class  of  variational  problems  originally  used  to  derive  bounds  for  the  effective 
parameters  of  composite  materials. 

In  Section  3  we  continue  the  examination  of  the  derivation  and  use  of  “bounds”  for 
the  effective  parameters  in  a  scattering  problem.  The  derivation  of  bounds  is  an  old 
subject  in  scattering  theory,  dating  from  early  in  this  century.  Nevertheless,  there  has 
been  some  very  interesting  recent  work  on  various  methods  for  deriving  bounds.  We 
discuss  this  work  briefly  in  Section  3,  as  it  applies  to  the  problem  of  scattering  from 
foliage.  As  with  homogenization,  the  class  of  bounds  discussed  applies  mostly  to  the 
low  frequency  case. 

The  work  on  the  use  of  “multiple-scale”  asymptotic  analysis  (homogenization)  for 
the  representation  of  scattering  processes  is  given  in  Section  4.  This  is  primarily  a  “low 
frequency”  theory,  valid  when  the  wave  length  of  the  incident  radiation  exceeds  the 
characteristic  dimensions  of  inter- particle  separation  in  the  scattering  medium.  The 
homogenization  method  does,  however,  provide  a  systematic  procedure  for  the  con¬ 
struction  of  high  order  approximations  to  the  “effective  parameter”  representations  of 
the  scattering  process.  These  approximations  may  be  evaluated  by  common  numerical 
procedures  when  the  scattering  medium  has  a  periodic  structure.  We  describe  this  case 
first,  before  turning  to  the  more  realistic,  but  less  tractable,  case  when  the  physical 
model  parameters  (permittivity,  conductivity,  and  permeability)  are  random  functions 
of  position  in  the  medium. 

In  Section  5  we  describe  a  procedure  combining  an  analytical  representation  for  the 
macroscopic  scattering  phenomena  from  a  “laver"  of  foliage  above  terrain  with  the  effec- 
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Figiirt’  1.1:  Angular  variation  of  tin*  i-yntliotif  pulso  response  from  a  sprurr  tree  at  S 
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tive  parameter  representations  of  the  microscopic  interactions  of  the  EM  radiation  and 
the  elementary  scatterers  (absorbing  centers)  in  the  foliage.  The  method  of  smoothing 
perturbation  is  used  to  obtain  a  deterministic  vector  wave  equation  for  the  average 
electric  field  scattered  by  the  random  “very  rough”  surface.  The  coherent  scattering  is 
characterized  by  reflection  coefficients  which  can  be  computed  easily.  The  technique  is 
applied  to  treat  surfaces  with  uniform  and  Gaussian  distributed  height  variations,  and 
the  results  are  compared  to  the  Fresnel  coefficients.  Only  small  differences  are  observed 
at  grazing  incidence;  however,  at  larger  incident  angles  substantial  differences  (20-50%) 
in  both  amplitude  and  phase  variations  are  obtained. 

In  Section  6  we  present  some  conclusions  and  recommendations  for  further  work. 


5 


2  T-matrix  Approximations  to  Multiple  Scattering 


In  this  section  we  use  the  transition  matrix  (T-matrix)  formalism  to  describe  multiple 
scattering  phenomena  in  a  region  containing  two  or  more  types  of  dielectric  scatterers 
arranged  in  a  random  geometry  in  a  stirrounding  background  medium.  Models  for 
scattering  from  a  medium  containing  a  two  (or  more)  classes  of  (randomly  oriented) 
scatterers  in  an  enveloping  medium  are  prescribed  in  terms  of  the  definition  and  deriva¬ 
tion  of  approximations  for  the  “effective  dielectric  constant”  resulting  in  a  family  of 
approximations  for  an  “effective  scattering  representation”  for  heterogenous  media. 


2.1  Basic  Framework 


Consider  a  region  O  of  space  occupied  by  a  homogeneous  (background)  dielectric 
medium  with  dielectric  constant  to  containing  elements  from  m  classes  of  scatterers 
with  dielectric  constants  ej,...,em.  Let  vtJ ,  t  --  1  —  1, 2, . . . ,  JV,,  be  the  sub¬ 

sets  of  O  occupied  by  scattering  elements  j  of  class  t.  The  dielectric  properties  of  be 
composite  material  in  O  are  described  by 

*71  N, 

f(r)  -  f0  4  -  fo)  \i'„(r)  (2.1, 

•  -1  J=i 


where 


\v,»  = 


1  r  e  r,j 


(2.2) 


(  0  r  £  r,j 

(See  Figures  2. 1,2. 2.)  Let  fi,  be  a  (dimensionless)  parameter  describing  a  characteristic 
dimension  of  elementary  scatterers  in  class  f  (e.g.,  the  radius  of  spherical  scatterers); 
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Figure  2.1:  Distribution  of  scatterers  and  plane  wave  incidence. 

and  let  p,  be  the  total  volume  fraction  of  O  occupied  by  scatterers  of  class  i  =  1, . . .  ,m. 
Suppose  a  constant  (in  amplitude  -  time  harmonic)  field  is  incident  on  the  region  0. 

We  wish  to  characterize  the  scattering  properties  of  the  composite  material  in  the 
limit  Ni  — *  oo,  f>,  — +  0,  with  p,,,t  =  constant.  We  are  particularly  interested 

in  evaluating  the  relative  roles  of  the  classes  of  scatterers  in  the  scattering  process,  e.g., 
in  terms  of  their  relative  densities  p,  in  the  cases  (a)  p,,  p,  — ♦  0  with  p,/p>  ~  0(1)  and 
(b)  Pi,P}  0  with  pi/pj  0. 

We  shall  sketch  the  development  of  a  family  of  effective  scattering  approximations 
for  the  material  using  the  multiple  scattering  framework  developed  originally  by  M. 
Lax  (1951,1952). 
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Figure  2.2:  Plane  wave  scattering  by  an  arbitrary  inhomogeneous,  anisotropic  object. 

2.2  Transition  Operator  Representations 

We  use  the  transition  operator  formalism  (Lax  1951,1952,1973)  as  used  by  (Lang 
1981)  and  (Kohler  and  Papanicolaou  1981)  among  many  others,  as  a  framework  for 
the  derivation  of  several  formulas  for  an  effective  dielectric  constant  for  the  composite 
medium.  Maxwell’s  equations  for  this  situation  may  be  written  in  the  form 

m  N, 

V  •  UoE)  +  £(e,  -  e0)  £  V  •  [x„, E]  =  0  (2.3) 

<  =  1  3- 1 

V  x  E  =  0 

with  the  boundary  condition 

<  E  >  —  E  hs  N,  — *  oo  (2 A) 
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where  <  •  >  is  expectation  and  E  is  the  constant  external  field.  The  constituent 
relationship 


D(x)  =  e(x)E{x) 


and  the  condition 


<  D(  )  >=  e’E 


defines  the  effective  dielectric  constant  for  the  composite  material. 


(2.5) 

(2.6) 


Rewriting  (2.3)  in  abstract  form 


(L0  +  M)E  =  0,  Vx£  =  0 


where 


Lo  —  V  •  (eo*) 

N, 

M  =  M\  *-•••  +  A/m, 

j=i 

Vn  -(€,•-  eo)-V(xu.,(-)] 


we  have 


E  +  L0'ME  =  F 

as  the  integral  form  of  (2.3).  F  is  chosen  so  that 

L0F  —  0,  <  E  >=E. 


(2.7) 


(2.8) 


(2.9) 


(2.10) 


The  transition  operator  (“matrix”)  T  is  defined  as  follows: 


L0'M  =  (Lo  +  M  -  M)~'M 
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=  [(Lo  +  M)(I  (L04 


(2.11) 


=  T(I  -  7T1 

with  T  =  {Lo  4-  M)_1M.  So  (2.9)  becomes 

[/  +  T(I  -  T)~1}E  -  F  (2.12) 

or 

E  =  (/  -  f)F 

The  condition  (2.10)  gives 

F  =  (/-  <  T  >)_1E 
E  =  (I  -  T)(I~  <  T  >)-JF 
Using  (2.5)(2.6)(2.15),  we  have 

c‘F  =  {<(>-<  (T  >)(/-  <  T  >)~XE 

as  the  definition  of  c'  in  terms  of  the  transition  operator. 

It  is  apparent  that  knowledge  of  T  completely  characterizes  the  scattering  properties 
of  the  composite  medium.  The  operator  T  is  related  to  the  scattering  amplitudes  of 
the  individual  particles,  and  to  their  interactions.  For  simple  (single)  dipole  scatterers, 
T  can  be  determined  from  the  polarizability  of  the  element;  and  so,  T  is  defined  in 
terms  of  quantities  of  physical  interest.  Other  properties  of  the  transition  operator  are 
described  in  (Lax  1951,1952,1973),  (Lang  1981),  and  (Kohler  and  Papanicolaou  1981). 
Our  analysis  of  the  representation  (2.15)(2.16)  is  based  on  (Kohler  and  Papanicolaou 
1981). 


(2.13) 

(2.14) 

(2.15) 

(2.16) 
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THE  TERM  INVOLVING  ONLY  Til) 

THE  TERM  INVOLVING  ONLY  T(2) 

THE  TERM  INVOLVING  T(2)CfT(l) 

THE  TERM  INVOLVING  T(l)0'T(2) 

THE  TERM  INVOLVING  T(l)O’T(2)0'T(l) 

Ficure  2.3:  Multiple  x-;» 1 1 *-ii n n  interpretation  of  the  T-matrix  formalism  for  two  scat- 
terrrs. 


I  I 


To  bring  out  the  interaction  of  the  various  classes  of  scatterers,  (see  Figure  2.2  we 
reconsider  (2.9) 

E  +  L^(M,+  •••  +  Mm)E  =  F  (2.17) 

which  may  be  written  as 


E  =  F  -  TjEi  -T2E2---TmEm  (2.18) 

with 

L0lM,  =  T,(I  -  T.r1 

Ti  =  (L0  +  Afi)-'Mi  (2.19) 

Et=(l-T,rlE 

since  E{  -  TtE,  =  E,  we  have 

rn 

Et  +  ^TkEk  =  F  (2.20) 

k?> 

Equations  (2.18 )-(2.20 )  express  the  field  in  terms  of  the  fields  and  transition  operators 
associated  with  each  class  of  scatterers  in  the  composite  medium. 


For  m  —  2  classes  of  scatterers  we  have 

T  =  (L0  +  M,  f  M2)-1  (Mi  +  A/2)  (2.21) 

=  {I  +  Tl2rlTl  +  (/  +  T2lr]T2 

where  the  pair  transition  operators  are 

Tt}  =  (Lo  +  M.V'Mj,  i  j  =  l,2,...,m  (2.22) 

This  representation  gives  the  formula 

c  ~[<c>  —  <  c(  /  +  7\2 )  1  Tj  >  -  <  ((I  +  T2\)  ^  7j  >  j  (2.23) 
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•[/-  <  (/  +  t12)-1t1  >  -  <  (/  +  T21)~'T2  >]-1 

for  the  effective  dielectric  constant  in  which  the  roles  of  each  class  of  scatterers  and 
their  interaction  is  explicit.  We  could  continue  this  process  to  make  the  role  of  each 
elementary  scatterer  of  class  i  explicit.  Note  that  (2.23)  contains  no  approximations. 


For  the  general  case  of  m  >  2  classes  of  embedded  scatterers  the  formula  is 


c  ~ 


;=i  j*i 


'-'£'E<v  +  T'1r'T.> 

.=1 Jjl 


-1 


(2.24) 


2.3  Average  T-Matrix  Approximation  (ATA) 

Consider  the  case  m  ~  2  for  a  moment.  Expanding  the  expression  for  the  T-operator, 
we  have 

T  =  Tx  +  T2  -  r12r,  -  TuT2  +  •  •  •  (2.25) 

neglecting  the  higher  order  terms  (TuTi, . . .),  and  substituting  in  (2.23),  we  obtain  the 
simple  approximation 

«■-!<*>-  E  <  eT<  >ll7  -E<Ti  >]-1  (2.26) 

i=i  ,=i 

called  the  average  T-matrix  approximation  (ATA)  (Lax  1951,1973)  (Kohler  and  Pa¬ 
panicolaou  1981).  The  expressions  <  Tj  >  and  <  >  in  (2.26)  diverge  in  typical 

cases  (Batchelor  1974)  (Gubernatis  1978),  so  the  ATA  must  be  used  with  care.  The 
usual  technique  is  to  employ  a  “cutoff”  argument,  limiting  attention  to  a  finite,  but 
large  region  of  space  and  ignoring  contributions  from  beyond  this  region. 
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It  is  possible  to  evaluate  (2.26)  in  the  limit  as  the  sizes  of  the  elementary  scatterers 
approaches  zero.  Suppose  each  class  consists  of  spherical  scatterers  having  radii  6,,i  = 
1,2  with 


Pi  - 


(2.27) 


the  volume  fraction  occupied  by  each  class.  Here  C{  is  the  average  number  of  sphere 
centers  per  unit  volume,  this  case  leads  to  the  approximation 


<  c  >  +Plfl 

[  fu-Q 

L  2<o  +<i  . 

1  +  p*>,  [g 

riSr) 

1  f  Px 

ii 

2f 

«ll  -C3 
,2<0  +*2  . 

(2.28) 


which  is  a  version  of  the  Claussius  -  Mossotti  formula  (Kohler  and  Papanicolaou  1981, 
p.  213)  (Landauer  1978).  In  general,  the  C'laussius  -  Mossotti  formula  is  viewed  as  a 
good  approximation  for  low  to  moderate  volume  fractions  pt. 


The  ATA  not  only  contains  divergent  terms  as  a  rule,  it  also  fails  to  account  for 
the  interactions  between  scatterers.  A  family  of  more  refined  approximations  has  been 
developed  to  include  interaction  effects  in  formulas  for  effective  scattering  representa¬ 
tions  which  are  nonetheless  computationally  feasible  (Lax  1951,1952,1973)  (Elliott  et 
al  1974). 


2.4  Coherent  Potential  Approximation  (CPA) 

Among  the  most  interesting  of  these  is  the  Coherent  Potential  Approximation  (CPA) 
which  essentially  involves  neglecting  the  difference  between  the  field  exciting  the  medium 
and  the  average  field.  Let  er  be  a  “reference”  dielectric  constant  to  be  selected,  and 
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consider  Maxwell’s  equations  rewritten  as 


V  l'rE  1 

«'  =  1  J=1 

m  N<  /  c  —  c  \ 

)'  E  =  0 

T  =  1  JTZ  1  1 


V  x  E  =  0 


We  may  write  this  abstractly  as 


with 


(Lr  +  Mr)E  =  0 


=  V-(cr-) 

2  A', 

<=) 


The  same  argument  as  before  yields  the  representation 

rn  N,  m  Nt 

*'  ~  l<  o-ZEoT’j  >1(1  -EE  <  is  >]- 

I  —  1  J- 1  1—1  J  — 1 


where 


^  =  l^r  +  Kj]_1  Kj 


The  CPA  is  based  on  choosing  (r  to  “optimize”  the  approximation  in  (2.32). 


Let  g{r)  be  a  smooth  function  and  consider  the  evaluation  of  T-- 

»’’V)  -  (Trg)(T) 


(2.29) 


(2.30) 


(2.31) 


(2.32) 


(2.33) 


(2.34) 
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so 


v . !«,»«)  +  V  ■  K<,  - 

=  V.|(., (^f)»l  (2-35) 

V  X  wij  =  0 

Let 

tr,N  =  \tr  +  --(fo  tr)]  =  fr  for  A,  large  (2.36) 

m  J\, 

and  w'J  =  w?  +  u>^  where 


(2.37) 

=  v  •  j(t,  -  #o).\f„y) 

V  -(t\N «•?)  +  V-[(t,  -eo)\v„«’ijl 

(2.38) 

-('irivr)*-' 

V  x  w\}  0  =  V  x  w'7J 


Let  {yj}  be  the  locations  of  the  centers  of  the  individual  scattering  elements.  For 
6,  small,  we  may  replace  g  on  the  right  in  (2.37)  by  y(yp)  =  9ij •  Using  this  and  the 
approximation 

('■N  *  cT  for  A,  large  (2.39) 


we  have 


to  —  t 


- — 


3tr  4-  t,  —  to  ' 
fo 


g u  •  (*  -  yj  )1 


lx  -  .Vui3 


1  >  *. 


(2.40) 


«’i  (-r) 


3fr  +  f,  tO 


9v 


lx  ;/,j  1  <  K 


16 


a  result  similar  to  (Kohler  and  Papanicolaou  1981,  sect.  8).  In  addition,  using  (2.39) 
in  (2.38) 

(.r  -  y)  ■  g(y ) 


where 


u'3V) 


<0  ~  t, 


1 


47rm  A’, 


1 


;v’/ 


l*  -  y I 
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-dy  4-  w3J 


3cr  4  f,  -  f0  m  A"; 

<0  "  f 


9ij  ■  (x  -  y,', ) 


k  -  y,j  I3 

l 


k  -  yH !  >  h 


nv(x  )  -  t  -  -  - .7<7u  ■ 

3fr  t  -  <o  niiV, 


k  -  y,T  I  <  <*< 


(2.41) 


(2.42) 


Combining  these,  we  have 


V(r,;?K.r)  - 

J  1 


•V, 


(r  -  y)  •  g(y) 


dy 


— f0-  -**  tfVfi  -  \,.„(.T))Vr 
3f r  t  ^ o  ^ 


r  -  y|* 

\9ij  ■  (*  -  y-4 ) 


_ Cl 

k  -  y/,  l3 


Ex,.,  k>9.,+  £“•>(*> 

3<r  4  *l  •  fo  ;V1  —  I 


(2.43) 


Arguing  as  in  (Kohler  and  Papanicolaou  1981,  sect.  8),  the  last,  term  on  the  right  in 
(2.43)  is  0(  Nt  1  )  on  the  average  as  Ar,  — *  oo.  We  take  the  expectation  and  choose  f, 

so  that 

(2.44) 


N, 


••  EEru 

,-i  j=i 


This  is  the  CPA.  This  requires 


fr  -V>.- 

3cr  ^  3fr  4  (,  c0 


(2.45) 


where  p,  =  4 jr^fc,  with  c,  the  average  number  of  scattering  centers  (of  class  i)  per  unit 


volume. 


17 


Using  this  choice  of  cr ,  the  effective  dielectric  constant  may  be  approximated  using 
the  CPA  by  evaluating  (2.32).  The  result  is 


*  =  €o  +  (  x  -  -  (1  -  p)e0  (\^1) 

-(!>•) 


(2.46) 


+  53  P'(' 


0=1 

<o  -  * 


3fr  4"  C)  —  fo 


where 


P  =  y^p, 


1  =  1 


This  result  should  be  contrasted  with  the  formula  from  Effective  Medium  Theory 
(EM A)  (Elliot  et  al  1974)  which  gives 

('  -  fo 


2f“  —  to 


'>+s>(£ 


o 


(2.47) 


The  most  accurate  formula  is  obtained  by  using  e‘  as  the  dielectric  constant  of  the 
“reference  medium.”  This  leads  to  the  Self-Consistent  (SCA)  formula 

<  6  >  +  ^->  =  1  P'f'  (  3r  • +e  *1»7,  ) 


(  — 


1  +  vm  ,  p  ( 

— ‘<-1  P'  \  3< '+«,-<«  / 


(2.48) 


The  results  from  several  approximations  (loss  versus  volume  fraction)  are  compared 
in  Figure  2.4.  It  is  clear  that  there  is  considerable  variation  among  the  methods. 


2.5  Selection  of  Trial  Fields  Based  on  Variational  Principles 

The  results  from  the  previous  subsection  can  be  put  in  perspective  by  considering  the 
problem  of  constructing  bounds  for  the  effective  parameters  representing  the  medium. 
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Imag  (EDC) 


EDC  (M= I )  VS  vol  free 


Figure  2.1:  ( 'ompnrison  <>f  the  effective  parameter  for  several  approximation  methods. 
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In  particular,  we  can  develop  the  process  of  selecting  a  “reference  background  medium” 
in  the  context  of  a  variational  problem. 

Let  £  be  a  curl  free  field,  and 

V  =  ^  <«£,£>,  V  x  £  =  0,  <£>=£  (2.49) 

where  V(-)  is  the  electrostatic  energy  in  £.  Then  minf  V(£)  is  achieved  by  E  satisfying 
Maxwell’s  Equations 

D(x )  =  f  {x)E(x) 

V  x  £  =  0,  V  D  =  0  (2.50) 

<  E  >=  E 

where  E  is  a  given  external  field.  Recall  <  D  >=  e"  E  defines  the  effective  dielectric 
constant.  If  £  is  any  field  with  V  x  f  =  0,  <£>=0,  then 

~  <  e£  •  £  >=  0  (2.51) 

Therefore, 

V{E)  =l-<tEE> 

=  ^  <  eE  •  E  >  +  ^  <  tE  •  (E  -  E)  > 

=  l-<iE>E  (2.52) 

=  ~  <  t  >  E  •  E 

2 

This  gives  the  general  upper  bound 

€'£•£<  <(£■£>  (2.53) 
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for  any  stationary  random  field  S  with  V  x  £  =  0,  <  €  >=  E.  Inequality  (2.53)  gives 
an  upper  bound  on  t  for  any  “trial  field”  £. 

To  obtain  a  lower  bound,  we  use  the  complementary  variational  principle.  Let 

\V(V)  =  l-  <  t~xV-V  >  (2.54) 

V  -V  =  0,  <  T>  >=  D 

The  minimum  is  again  achieved  by  Maxwell’s  equations.  This  provides  the  bound 

(e')~lDD<  <t~xV-V>  (2.55) 

for  any  trial  field  T>  with  V  •  T>  —  0,  <.V  >~  D. 

To  construct  suitable  trial  fields,  we  use  the  transition  operator  -  multiple  scatter¬ 
ing  formalism.  Recall  the  formulation  of  Maxwell’s  equations  for  the  multicomponent 
medium 

m  N, 

V  ■  (e0E)  +  £  £  V  ■  [(«,  -  to)\v>JE]  =  0  (2.56) 

J= i 

V  x  E  =  0 

which  we  write  in  abstract  form  as 

(Lo  +  M)E  =  0,  V  x  E  =  0  (2.57) 

As  before 

E  +  Lo1ME  =  F  (2.58) 

where  F  satisfies  LoF  =  0,  <  E  >—  E.  In  terms  of  the  transition  operator  T  = 
(Lo  +  M)~x M  we  may  write  this  as 

E  =  (I-  T)F 
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(2.59) 


=  V  Er')F 

t=i 

m  N, 

!=1  j=l 

Based  on  this  we  define  the  trial  field 

S  =  F  El tT;F  (2.60) 

■=i  j=i 

with  F  a  constant  to  be  chosen.  Writing  out  the  operator  T,j  explicitly  we  have 

+  X„.»F}  (2.61) 

This  trial  field  is  curl  free.  We  choose  F  so  that  <  £  >—  E. 


|(.t  -  ytj)  •  F 

I*  -  y.j  I3 


We  shall  continue  this  analysis  in  the  next,  section  which  discusses  general  procedures 
for  constructing  bounds. 
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3  Bounds  for  Effective  Parameters 


In  the  previous  section  and  in  the  next  we  present  methods  for  representation  and  com¬ 
putation  of  effective  parameter  models  for  heterogeneous  media  in  terms  of  statistical 
characterizations  of  the  media.  These  methods  lead  to  representations  for  the  effective 
parameters,  e.g.,  the  effective  dielectric  constant  and  conductivity,  which  may  be  diffi¬ 
cult  to  evaluate.  In  this  section  we  continue  the  examination  of  the  derivation  of  bounds 
on  the  effective  parameters  as  a  means  of  quickly  providing  good  approximations  to 
the  behavior  of  the  medium.  We  elaborate  on  the  ideas  introduced  at  the  end  of  the 
previous  section  for  construction  of  bounds  based  on  variational  procedures.  We  also 
consider  the  construction  of  bounds  based  on  the  method  of  analytical  continuation 
introduced  by  Bergmann  and  extended  by  Willis,  Golden  and  Papanicolaou. 

Approximations  to  the  effective  permittivity  (or  conductivity  in  thermal  systems) 
have  been  derived  by  three  basic  methods.  The  original  work  of  Haskin  and  Shtrik- 
man  (1962)  was  based  on  variational  principles,  and  these  have  been  reformulated  in 
recent  years  to  provide  one  of  the  basic  methods  for  deriving  aproximations  (Milton 
1981a, b).  The  second  method  for  deriving  bounds  is  the  method  of  “compensated  com¬ 
pactness”  (Tartar  1986).  The  third  method  involves  the  use  of  complex  analysis  and 
“representation  formulas”  (Bergmann  1978,  Golden  1986,  Golden  and  Papanicolaou 
1983). 

There  has  been  a  great  deal  of  work  in  this  area,  not  only  in  electromagnetics, 
but  also  in  conductivity  and  elasticity.  This  work  dates  from  the  early  years  of  this 
century  (Bergmann  1978),  and  it  includes  some  very  recent  work  (e.g.,  (Ericksen  1986)). 
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We  shall  describe  the  applications  of  this  work  for  the  description  of  scattering  and 
absorption  of  EM  radiation  by  foliage  covered  terrain. 


3.1  Definition  of  the  Effective  Parameters 

The  starting  point  is  a  definition  of  the  effective  parameter,  e.g.,  the  effective  permit¬ 
tivity.  We  shall  use  a  slightly  more  general  model  than  the  one  used  in  the  previous 
section.  Suppose 

f(x,u>)  =  +  •••  +  £m.\m(u>)  (3.1) 

is  the  dielectric  “constant”  of  a  composite  medium,  where  €,,  t  =  1, . . .  ,m  is  the  (com¬ 
plex)  dielectric  constant  of  medium  i,  and  the  indicator  function  is  1  for  all  samples 
u ;  €  fi  which  have  medium  i  at  point  x,  and  zero  otherwise.1  We  use  {Q.,Jr,V)  to 
denote  the  probability  space  on  which  the  st  atistics  of  the  phenomena  are  defined.  The 
quantity 

Pi=  /  (3-2) 

J  n 

is  the  volume  fraction  occupied  by  the  material  of  type  j  =  1,2, ...  ,m. 

The  spatially  dependent  random  conductivity  and  magnetic  permeability  are  sim¬ 
ilarly  defined.  In  our  case  these  functions  depend  on  the  random  geometry  of  the 
scattering  foliage  which  may  be  available  only  in  terms  of  second  order  statistics.  The 
bounds  obtained  for  the  “effective  parameter”  approximations  of  these  functions  should 
not  depend  on  statistics  beyond  this  level. 

‘Thus,  we  are  assuming  m  classes  of  scatterers  as  in  the  previous  section.  We  shall  not  have  to  deal 
with  the  individual  scattering  elements  at  this  point. 
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If  E(x  ,u>)  and  D(x, u>)  are  the  electric  and  displacement  field  vectors  (stationary  in 
time),  then 

Dt(x,u>)  =  tE,(x,u>)  (3.3) 

V  ■  D{x,v)  =  0 

V  x  E(x,u>)  =  0 

and 

f  P(dw)E(x,u>)  =  E 
J  o 

provides  the  boundary  condition,  where  E  is  a  constant  field  incident  on  the  region.  By 
normalizing  the  field  magnitudes  and  redefining  the  coordinates,  we  can  take  E  =  e*, 
the  kth  unit  vector  in  9?3.  Let  Ek  and  Dk  denote  the  corresponding  fields.  Then  (3.3) 
may  be  rewritten 

d 

Dk(x,u>)  =  £  «,>(*, u>)£*(*,u>)  (3.4) 

0=i 

The  effective  dielectric  constant  may  be  defined  as 

=  [  P(du)Dk(u})  (3.5) 

Jo 

In  (Golden  and  Papanicolaou  1983)  it  was  shown  that  this  essemble  average  coincides 
with  the  usual  definition  involving  a  volume  average. 

Equation  (3.5)  may  also  be  written  in  the  symmetric  form 

=  (3.6) 

where  the  bar  denotes  complex  conjugate.  This  expression  suggests  the  interpretation 
of  f'  as  the  dielectric  constant  of  a  fictitious  homogeneous  medium  (for  the  volume 
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V)  which  provides  the  same  value  for  the  electrostatic  energy  stored  in  the  volume 
occupied  by  the  heterogeneous  medium.2 


3.2  Variational  Methods  for  Computing  Bounds 

In  the  last  section  we  introduced  the  variational  approach  to  the  computation  of  bounds 
on  the  effective  dielectric  constant.  In  this  subsection  we  shall  develop  this  approach 
in  detail  -  providing  an  assessment  of  the  utility  of  the  method  for  the  treatment  of 
scattering  and  absorption  phenomonena  in  foliage  and  several  formulas  for  effective 
parametric  representations. 

The  starting  point  is  the  pair  of  variational  principles: 


f "  =  min  <(:£■£> 

Vxr=n 

<£>~t 

(3.7) 

(e')_1  =  min  <-T>-V> 

v  r>=s<i  £ 

<Vy~t 

(3.8) 

Here  e  is  a  unit  vector  (normalized  external  field),  and  £,V  run  over  stationary,  square 
integrable  random  fields.  The  bracket ts  <  •  >  as  usual  denote  ensemble  average. 

The  minima  are  attained  by  the  solutions  of  Maxwell’s  equations  (3.3).  Note  also 
that 

e*  =<  T>  •  e  >  =  <(£■  e  >  =  <  e£  •  £  >  (3.9) 

2The  volume  average  form  is  e”t|£o|3  =  [/v  e(r)|£(r)|::dV']  /ineas(V)  where  V  is  the  representative 
region,  meas(V)  is  its  volume,  and  E»  —  (/^  E(r)dV\  /rneas(V)  is  the  average  (incident)  field. 
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because  <  eS  •  E  >=  0  for  any  test  field  E  with  <  E  >=  0  and  V  x  E  =  0.  It  is  also 
true  that 

<  ~(t-yxv-(e')-xv>=  (emr2  <x-v-v>  (3.10) 

c  e 

-  (O-2  <(EE> 

—  ( f * )  2  <  T>  •  e  >—  (c")_l 
Hence  the  two  variational  problems  are  consistent. 

The  construction  of  bounds  for  e~  reduces  to  selecting  suitable  trial  fields  E,T>  to 
insert  into  the  variational  problems.  Choosing  E  =  e,T>  =  e  leads  to  the  Wiener  bounds 

(<  ~  >)_1  <  e"  <  <  (  > 

( 

known  since  early  in  the  century  (Bergmann  1978). 

As  pointed  out  by  Kohler  and  Papanicolaou  (1982),  there  are  two  basic  considera¬ 
tions  to  the  construction  of  bounds: 

1.  Utilize  the  underlying  structure  of  the  medium  (“its  geometry”)  in  the  analysis, 
e.g.,  dielectric  spheres  embedded  in  a  homogeneous  medium,  to  produce  geometry 
dependent  bounds ;  or 

2.  Produce  bounds  which  are  geometry  independent  for  materials  in  which  the 
components  have  a  complex  infrastructure. 

The  Claussius-Mossotti  formulas  are  examples  of  bounds  of  the  first  kind.  The 
results  obtained  using  the  Coherent  Potential  Approximation  in  the  previous  section 
are  examples  of  geometry  dependent  bounds. 
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3.2.1  Geometry  Independent  Bounds 

As  in  the  previous  section,  we  introduce  a  parameter  fr  which  is  the  dielectric  constant 
of  a  “reference  medium”  and  rewrite  Maxwell’s  equations  as 

V-(«r£)  + V  •((<-€,)•£]  =0,  Vx£  =  0  (3.11) 

Introducing  the  abstract  operator 

r£  =  -V(-A)_1(V  •  E)  (3.12) 

we  can  rewrite  (3.11)  as 

E  +  T(r1(e  -  tr)E-  <  r‘(«  -  fr )E  >)  =  E  (3.13) 

(which  actually  encompasses  all  of  Maxwell’s  equations).  As  pointed  out  in  (Kohler 
and  Papanicoloau  1982),  the  operator  T  is  a  projection  on  the  family  of  zero  mean, 
stationary  random  fields  into  the  subspace  of  curl-free  random  fields.  Formally, 

r2  -  r 

In  addition,  if  £  is  a  stationary  random  field  which  has  zero  mean  and  is  isotropic ,  i.e., 
the  correlation  function  satisfies 


<  £,(: r  +  y)£j(y)  > =  H,,(|x|) 


then 


<T££>  = 


1 

3 


<  £  ■£  > 


(3.14) 
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To  associate  (3.13)  with  a  variational  problem,  we  follow  (Kohler  and  Papanicolaou 
1982)  and  introduce  the  polarization  field 

P  =  (;\e-er)E  (3.15) 

which  satisfies 

er(f  -  er)-1P  +  r(P-  <  P  >)  =  E  (3.16) 

Equation  (3.16)  is  the  Euler  equation  for  the  quadratic  functional 

V(V)  =■-  '  <  er(e  -tT)-xV-V  >  (3.17) 

+  \  <  rep-  <  V  >)-(V-  <V  >)>  -  <e-V  > 

Li 

When  er  <  e(x,u)),\/(x,u)),  then  U{V)  achieves  its  minimum  at  P  the  solution  of  (3.16) 
and 

IHP)  =  1  (i  -  9  (3.18) 


Thus, 


-  ^1  -  —  j  <  U(V),  when  (e  -  er)  >  0 
for  any  trial  field  V.  This  bound  can  be  rewritten  as 


cr(l  -  2 U(V))  <  when  (e  -  er)  >  0  (3.19) 


When  («  —  €r)  <  0  the  quadratic  form  U{V)  is  negative  definite.  This  provides  the 
upper  bound 


e*  <  «r(l  -  2U(V)),  when  (t  -  er)  <  0 


(3.20) 


The  objective  of  the  analysis  is  to  choose  trial  fields  V  to  make  the  bounds  ( 3.1 9)(  3.20) 
as  tight  as  possible.  Rewriting  (3.15)  as 

p  -  LZJl  jf  _  r(P-  <  P  >)j  (3.21) 

there  are  several  different  ways  to  proceed. 


First,  if  the  medium  is  nearly  uniform:  |£t|  =  |e  —  er\  is  small  for  some  er,  then  we 

can  expand  P  in  a  power  series  in  Se,  and  produce  trial  fields  by  truncating  the  series; 

_  Se  Se  I  Se  Se  \  , 

V  =  —  e - 1'  -e-  <  —  e  >  (3.22) 

<r  t r  V  (r  (r  ) 

This  case  has  been  treated  extensively;  see  (Bergmann  1978)  (Milton  1981a, b)  (Golden 
and  Papanicolaou  1983). 


Second,  if  Se  is  not  small,  but  the  medium  has  a  host-inclusion  geometry.  Then 
the  multiple  scattering  formalism  can  be  used  to  derive  bounds.  This  was  illustrated 
in  the  previous  section.  We  shall  consider  it  further  in  the  next  subsection. 

In  the  remainder  of  this  subsection  we  shall  consider  the  third  alternative  -  Se  is 
not  small  for  any  cr  and  the  medium  is  not  of  the  simple  host-inclusion  type.  However, 
the  medium  (the  random  function  e(:r,u>))  is  statistically  isotropic. 


Introducing 


To  -  -  /and  r,  =  F  -  To 


we  have 


<  r,£,  •  £"2  >—  0.  <'  rfj  •  £*2  >  —  ~  <  £]  •  £2  y 
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for  any  fields  £, ,  i  =  1,2  that  are  isotropic.  Using  this  notation,  we  can  rewrite  (3.15) 
in  the  form 


p-  <p>  +-r0(P-  <  p  >)  +  — r,(P-  <  p  >)  =  —e—  <  p  > 

( r  fr  (r 


or 


r*\ 

1 

1% 

_ i 

(  St  \  1 

St 

/+  J  +  -r0  — rt 

( p -  <  p>)=  /  +  -r0 

— e-  <  P  > 

L  '  tT)tT 

V  £r  / 

(3.23) 


Defining 


+  77r' 


we  can  write  (3.23)  in  the  form 


p-<p>=P~  {i+7,r°y'(z'-<p>) 


(3.24) 


We  can  choose  trial  fields  by  taking  successive  terms  in  this  series.  The  simplest  is 

(3.25) 


v-  <  v  >=  ^/  +  {~e~<v>) 


Taking  averages  in  (3.25),  and  using  the  definition  of  T0,  we  have 


<  V  >  = 


4*fef ']"<("£)"! 


>  e 


and 


U(V)=  ~-<V>-e 


(3.26) 


(3.27) 


Using  (3.26)(3.27)  in  (3.19)(3.20)  we  have 


!■* 

/ 

St  \ 

-1 

I ,  St) 

1  St  ) 

H 

X  +  ZZ) 

> 

<  V  +  iZ) 

A 

1  ^ 

<  (’ 


(3.28) 
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when  8(  >  0 


and 


e"  <  (r  (  1  + 


<  1  + 


3  ej 


-l 


1  -1 


> 


<  1  + 


8e 

3^ 


1  6( 

—  > 
(r 


(3.29) 


when  St  <  0 


Note  that  these  bounds  apply  to  general  media  which  have  a  statistically  isotropic  struc¬ 
ture. 


When  the  medium  has  only  two  components  -  t(x,u ;)  takes  values  e0»fo  with  0  < 
e0  <  ei  <  oo,  then  inequality  (3.28)  is  optimized  by  taking  cr  =  eo‘.  and  inequality  (3.29) 
is  optimized  by  taking  cr  =  «j.  This  leads  to  the  classical  bounds  of  Hashin  and 
Shtrikman  (1962) 


«o  + 


Pi 


<  ('  <  ft  + 


P  o 


(l/(<i  —  Co)  +  ^o/3to  (l/(fo  —  (i )  4-  pi/3ci 

As  before,  pi,i  =  0,1  is  the  volume  fraction  occupied  by  material  of  type  t. 


(3.30) 


By  takng  additional  terms  in  the  expansion  (3.24)  one  can  obtain  better  trial  fields. 
Note,  however,  that  each  succeeding  trial  field  requires  even  more  statistical  informa¬ 
tion.  For  example,  the  next  trial  field  corresponding  to  keeping  terms  n  =  0, 1  in  (3.24) 
requires  second  order  statistics  (e.g.,  correlation  functions)  on  the  random  media. 


3.2.2  Geometry  Dependent  Bounds 

Now  suppose  that  material  is  of  the  host -inclusion  type  with  a  specific  geometry, 
spheres  with  dielectric  constants  e =  1 . m  embedded  in  a  uniform  background 
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material.  Thus,  let  y,j,i  =  1  =  1 , . . . ,  TV*  be  the  (random)  locations  of  centers 

of  the  spheres,  each  of  radius  6.  Assume  that  the  spheres  do  not  overlap.  Then 


e{z,u>) 


|t  -  y,j|  <  6,  for  some  i,j 
I*  -  y.j|  >  6,  for  all  i,  j 


(3.31) 


As  before,  px  —  ^7r£3c,  is  the  volume  fraction  occupied  by  the  spheres  or  class  i.  where 
c,  is  the  average  number  of  sphere  centers  per  unit  volume. 


To  construct  suitable  trial  fields,  we  use  the  transition  operator  -  multiple  scatter¬ 
ing  formalism.  Recall  the  formulation  of  Maxwell’s  equations  for  the  multicomponent 
medium 

m  N, 

v  •  (to E)  +  EE  V  •  [(<.  -  «, )x.„E]  =  0  (3.32) 

.  =  1  j=l 

V  X  E  =  0 

which  we  write  in  abstract  form  as 

{L0  +  M  )E  =  0,  Vx£  =  0  (3.33) 

As  before 

E  +  L~'ME  =  F  (3.34) 

where  F  satsifies  L0F  —  0,  <  E  >—  E.  In  terms  of  the  transition  operator  T  — 
(L0  ■f  M)  lM  we  may  write  this  as 


E  -  (I  -  T)F 

=  (/-Er,)F  (3.35) 

1-1 
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m  Ar, 

i=l 3  =  1 


Based  on  this  we  define  the  trial  field 


N, 


f  =  (3-36) 

«=i  j=i 

with  F  a  constant  to  be  chosen.  Writing  out  the  operator  TX}  explicitly  we  have 

*  r  r(*-y,,)-Fl 


+  XvtJ(*)F  >  (3.37) 


I*  -  vav 

This  trial  field  is  curl  free.  We  choose  F  so  that  <  £  >=  E.  If  we  take  {y,j }  to  be 
Poisson  process  with  parameter  c,  and  if  we  assume  a  suitable  cutoff  region  so  the  series 

$"f3cF  (3.38) 


converges,3  then 

m  / 

1=1  x 

Hence, 

r=(  i-E 

'  i=i 

Now  let 

€.  -  fO 

7.  =  ^ - 

2f,  -  e0 

-1 


(3.39) 


,  t  =  1 , . . . ,  m 


.r  •  E 


3(2) 


then  (3.37) 

/  m  > 

£  =  (i  -  Y  p'~)' 


- 1 


i=i 


m  A’, 


E  -  Y  Y{EU](X  ~  V'J  •( 1  -  \^x  ~  .Vo')) 


(3.40) 


.=i  j= i 


+  E\7\t  -  ytj)\s{ x  -  )}] 


3This  is  the  case  of  the  ATA  divergence. 
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where  is  the  characteristic  function  of  the  sphere  of  radius  6  centered  at  the  origin. 


For  the  case  m  =  1  (one  class  of  embedded  scatterers),  this  formula  yields  the  upper 
bound  (Kohler  Papanicolaou  1981)  (using  the  Poisson  distribution  for  the  locations  of 
the  scattering  centers) 


.  1  +  up  -  2-f(p  4-  /tp  4-  *p2)  +  72[5p  +  7 up2  4 -  Kp  +  p2  +  Kpz ] 

(1  -  7p)2 

where  7  =  7 \,k  =  Ci/c0.k  =  k  -  1,  and  p  =  px. 


(3-41) 


A  lower  bound  can  be  obtained  by  a  similar  procedure,  starting  from  the  com¬ 
plementary  variational  principle.  The  trial  field  is  computed  again  by  expanding  the 
expression  for  7^ .  The  resulting  field  is 

®  =  O  ~  1 1  {(l  -  (771^)  -  27«.,(x)g|  (3.42) 

which  satisfies  <  V  >—  D,  assuming  the  elementary  scatterers  do  not  overlap.  To  have 
<  V  >=  D,  we  must  have 

g  =  {i  +  2]T7,p,rii} 

i  =  l 

In  the  case  m  =  1,  direct  computation  of  the  average  <  eV  ■  V  >  yields  the  estimate 


(O  < 


1  1  -f  K^p  -  27(p  +  k ip  -  2kip2)  +  7*(8p  +  16«ip2  +  4 p2  4-  4K1p3) 


to 


(1  4-  27 p)2 


(3.43) 


where  7  is  as  above,  and  =  k  1  -  l,fc  =  Ci/cq. 


Improved  bounds  may  be  obtained  using  the  Coherent  Potential  Approximation  to 
compute  the  expression  for  Ttj .  It  is  not  necessary  to  repeat  the  computation.  Rather 
this  may  be  accomplished  simply  by  introducing  the  reference  dielectric  constant  er. 
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The  trial  fields  are  constructed  by  assuming  that  the  exterior  dielectric  constant  is  er 
and  that  the  dielectric  constant  inside  the  scatterers  is  er  +  e,  -  c0  rather  than  c,.  This 
changes  the  trial  fields,  altering  the  parameter  to 


Now  7 i  varies  on  the  interval 


7l-  - 


t,  -  to 


3cr  4-  e,  --  (q 


ki  —  1  ki  —  1 

4k,  -  1  <  7’  <  it,-  +  2 


where  fcj  =  ti/to  and  we  assume  0  <  e0  <  e,. 


(3.44) 


(3.45) 


3.3  Bounds  Based  on  the  Method  of  Analytic  Continuation 

In  this  subsection  we  shall  examine  the  analytical  continuation  approach  to  deriving 
bounds,  indicate  its  range  of  applicability,  and  discuss  techniques  for  adapting  the 
method  to  treat  scattering  and  absorption  from  foliage.  The  method  is  based  on  writing 
the  effective  dielectric  const  am  as  a  function  of  the  (complex)  dielectric  constants  of 
the  underlying  component  media.  Hence,  it  is  represented  as  an  analytic  function  of 
a  complex  variable,  and  the  properties  of  such  functions  can  be  used  to  derive  bounds 
on  the  effective  parameter. 

3.3.1  Bounds  for  Two  and  Three  Component  Media 


If  we  substitute  (3.1)  into  (3.5)  and  divide  by  then  we  obtain 

rn;fc(hi,...,/rjv-i)  =  =  Ju  ^(rfu;'  ^  II  +  X/vj  Ef(u> 


(3.46) 
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Here  hi  =  The  key  idea  in  the  work  of  Bergman  (Bergmann  1978)  is  to  regard 

m(-)  as  a  function  of  the  complex  variables  hi, . . .  ,/i/v-i  in  CN~l.  This  permits  one  to 
use  (spectral)  representation  formulas  for  m(-)  which  allow  the  derivation  of  bounds 
based  on  the  extreme  points  of  certain  sets  of  measures.  By  expanding  the  represen¬ 
tation  formulas  about  the  case  of  a  “homogeneous”  medium  hj  —  1  ,j  =  —  1, 

it  is  possible  to  obtain  a  characterization  of  the  underlying  measures  in  terms  of  their 
“moments.”  These  moments  are  functions  of  increasing  amounts  of  information  on 
the  statistical  properties  of  the  medium.  For  example,  if  only  the  volume  fractions  pj 
from  (3.2)  are  known,  then  only  the  first  term  in  the  representation  can  be  computed 
explicitly.  This  leads  to  the  classical  Wiener  bounds: 


For  two-component  media  - 

(—  +-  <  *  <  Piti  4-  p2«2 

Ve,  c?  / 


(3.47) 


For  three-component  media  - 

(—  4-~  +  <  t  <  PiCj  +  p2f2  +  P3C3 

v  fl  (2  f3  ' 

These  bounds  are  achieved  by  parallel  plane  configurations  of  the  materials. 


(3.48) 


If  additional  information  is  available  on  the  statistical  description  of  the  medium, 
then  more  sophisticated  bounds  can  be  obtained.  For  example,  for  a  two-component 
medium  which  is  statistically  isotropic,  the  first  moment  of  the  spectral  measure  can 
be  computed  (in  the  expansion).  This  leads  to  the  Hashin-Shtrikman  bounds  (assume 
without  loss  of  generality  €]  <  f2) 

Pi  ^  ^  ,  Pi 


fi  + 


l/(«2  --  Cl  )  4  Pt/3ci 


S  f  5;  f2  T 


l/(ei  ~  €2)  4-  p2 /3e 2 


(3.49) 
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o.o  0.1  0.2  0.3  0-4 


Vol  froc  of  scol  terer  type  one 


Figure  3.1:  Effective  permittivity  from  EM  A  versus  bound  estimates. 

In  Figure  3.3.1  we  show  a  typical  result  illustrating  the  relationship  between  the 
upper  and  lower  bounds  and  the  effective  permittivity  (real)  versus  volume  fraction  as 
computed  by  the  T-matrix  EM  A  method. 
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3.3.2  Bounds  for  Multicomponent  Media 


There  is  a  substantial  difficulty  in  passing  from  two-component  materials  to  ^-component 
materials  with  N  >  3  when  using  the  analytical  continuation  method.  Basically,  one 
must  deal  with  functions  of  more  than  one  complex  variable  (/»i, .  • .).  The  precise 
difficulty  arises  from  the  fact  that  the  extreme  points  of  the  associated  set  of  spectral 
measures  are  not  known.  Prior  to  the  paper  (Golden  Papanicolaou  1985)  there  was 
no  systematic  method  for  treatment  of  multi-component  media  (using  the  continuation 
method).  In  the  case  of  three-component  media,  it  is  possible  to  circumvent  this  prob¬ 
lem  by  imposing  a  linear  relationship  among  the  complex  parameters,  and  expanding 
the  spectral  measure  about  the  (linear)  parameter.  In  (Golden  Papanicolaou  1985)  a 
technical  hypothesis  was  introduced  which  permitted  recovery  of  the  Hashin-Shtrikman 
bounds  for  three-component  media.  Assuming  Ci  <  c2  <  e 3  we  have 


ei 


*  -  S’"  th;  *  i) 


3c3 


(3.50) 


These  bounds  are  “optimal”  for  materials  consisting  a  mixture  of  spheres  of  all  different 
sizes  of  ci  and  e2  materials  each  coated  with  c3  material  in  the^ppropriate  volume 
fraction.  Additional  bounds  were  obtained  in  (Golden  1986)  for  complex  permittivities 
(thus  including  the  absorption  of  energy  in  the  material). 


The  analytical  continuation  method  is  a  powerful  technique  for  derivation  of  bounds 
on  the  material  properties  of  multi-component  media;  and  it  is  useful  to  examine  its 
potential  for  the  evaluation  of  scattering/absorption  phenomena  in  the  interaction  of 
EM  radiation  and  foliage.  The  following  points  are  important: 
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•  The  method  is  based  on  the  assumption  that  the  medium  is  nearly  homogeneous; 
that  is,  the  normalized  permittivities  t j/tN  ~  1.  This  is  a  restrictive  assumption 
for  foliage,  and  it  may  limit  application  of  the  approximation  to  certain  kinds  of 
foliage  (mixtures  of  grasses,  etc.). 

•  While  the  case  of  two  component  media  is  reasonably  well  developed;  the  case  of 
three-component  media  is  less  clear.  This  is  the  key  case  if  one  is  to  understand 
the  interaction  of  scattering  effects  in  an  environment  of  two  different  types  of 
scatterers  in  a  uniform  background.  The  limitations  on  the  method  are  both 
technical  and  physical. 

—  The  technical  limitations  are  mainly  due  to  the  incomplete  characterization 
of  the  set  of  spectral  measures  for  functions  of  two  complex  variables. 

-  The  physical  limitations  are  more  fundamental.  The  method  does  not  ac¬ 
count  for  the  geometry  of  multiple  scattering  processes.  That  is,  while  it 
is  capable,  in  principle,  of  including  higher  order  statistical  functions  in  the 
expansion  for  the  effective  permittivity,  the  method  develops  the  bounds  in 
a  fashion  which  is  well  removed  from  the  physical  properties  of  the  medium. 
The  extreme  points  of  the  set  of  measures  which  are  used  to  determine  the 
bounds  are  singular  measures  (“delta  functions")  which  may  not  correspond 
well  to  the  distributions  of  scatterers  in  typical  foliage.  It  would  be  better  to 
develop  the  bounds  in  a  (restricted)  set  of  measures  which  are  representative 
of  foliage. 
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4  Homogenization  and  Multiple  Scattering 


In  this  section  we  turn  our  attention  to  an  alternative  method  for  representation  of 
multiple  scattering  effects  in  heterogeneous  media.  The  method  is  a  variation  of  the 
“homogenization”  procedure  which  has  been  used  widely  in  mathematical  physics  and 
engineering  to  develop  effective  media  approximations.  This  method  has  promise  for 
for  several  reasons: 

1.  It  leads  more  naturally  than  the  T-operator  formalism  to  a  compact  represen¬ 
tation  for  an  effective  parameter  (complex  dielectric  constant)  representation  for 
the  effective  medium  approximation. 

2.  The  resulting  representation  for  the  effective  dielectric  constant  includes  the  in¬ 
teraction  of  microscopic  effects  (multiple  scattering)  explicitly. 

3.  The  underlying  analysis  applies  (however,  with  significant  differences  in  detail) 
to  both  periodic  and  random  media.  In  the  former  case  it  is  possible  to  solve  the 
equations  for  the  approximation;  in  the  latter  it  is  necessary  to  develop  approxi¬ 
mations  which  use  physically  measurable  quantities  (second  order  statistics  and 
correlation  functions). 

4.  One  can  prove  convergence  of  the  scaled  model  to  the  “homogenized”  model  in 
both  the  periodic  and  random  case. 

5.  The  representation  provides  a  (formal)  basis  for  the  systematic  construction  of 
a  sequential  approximation  to  the  effective  dielectric  constant,  including  “higher 


41 


order”  expansions  for  the  effective  parameters  in  terms  of  the  small  parameter, 
at  least  in  the  periodic  case. 

We  wish  to  undertake  a  systematic  investigation  of  such  approximations  in  the 
context  of  scattering  from  foliage  covered  terrain,  and  the  multiscale  (homogenization) 
method  offers  a  more  general  setting  for  such  a  comparison  than  does  the  T-matrix 
formalism,  at  least  in  the  frequency  range  to  which  it  applies. 


4.1  A  General  Model 


To  illustrate  the  ideas,  consider  the  following  general  model:  Let  O  C  be  a  region  in 
which  t(x),/j(x)  and  <r(x)  are,  respectively,  the  dielectric  tensor,  magnetic  permeability 
tensor,  and  conductivity  tensor  (3x3  matrices)  of  the  material  in  region  O.  In  O  the 
(vector- valued)  electric  and  magnetic  fields  satisfy 


d 

e(a:) 

0 

’  E(x)  ' 

di 

0 

//(a-) 

H(t) 

E(1,x) 

H{t,x) 


(4.1) 


+ 


<r(x)  0 

0  0 


E(t ,*) 

H(i,x) 


E( 0,r)  =  Eq( x),  H(Q,x)  -  H0(x) 


V  ■  [cJ^ol  =  O.V>ff0]  -  0 


(4.2) 

(4.3) 


where  V x  is  the  curl  operator  and  V-  -  div.  Note  that  (4.3)  implies 


V  •  [c£'(/,.t)J  =  0,V  ■  !/(//(t.r)j  =  0.V7  >  0 


■12 


It  is  necessary  to  assume  that  e(a-),^(x)  and  cr(x)  are  symmetric ,  and  that  f(x),/i(x) 
are  positive  definite  matrices. 


using 


V  x  E  = 


0  -  —  — 
dx 3  dx2 


d 

dxi 

—  —  0 

dx2  dxi 


*  0 
a 


’  E1  ' 

e7 

E 3 

(4.4) 


we  see  that  (4.1)  may  be  written  in  the  abstract  form 


where 


y(<,.r)  = 


a*, 

E(/,x) 

#(/,x) 


(4.5) 


G  ft6 


(4.6) 


>t°(x)  = 


£(*)  0 

0  fi(x) 

the  Ap(t,i)  are  skew-symmetric,  6x6  matrices  with  0’s  and  l’s  as  c  clients.  and 


B(x)  = 


cr{x)  0 

0  0 


G  »' 


6X6 


The  symmetry  assumption  means  that  (4.5)  is  a  symmetric  •  hyperbolic  system  of 
a  type  common  in  mathematical  physics.  The  (asymptotic)  analysis  of  systems  of  this 
type  with  periodic  coefficients  depending  on  a  small  (dimensionless)  parameter,  e.g., 

4”(x/i»),  Al[x/v),Bv{x/v) 
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(for  0  <  v  <<  1)  was  carried  out  in  (Bensoussan,  Lions,  and  Papanicolaou  1978)  and 
in  (Sanchez-Palencia  1980)  using  multi-scale  methods  (homogenization).  The  periodic 
case  represents  media  with  regularly  imbedded  components. 

The  “inhomogeneities”  in  the  media  which  govern  the  scattering  of  EM  radiation 
from  foliage  covered  terrain  are  randomly  distributed,  and  so,  the  homogenization 
(multi-scale)  method  must  be  extended  to  treat  this  case.  A  general  theory  for  partial 
differential  equations  with  random  coefficients  was  developed  in  (Papanicolaou  and 
Varadhan  1979).  This  work  forms  the  basis  for  our  treatment  of  the  random  scattering 
(and  absorption)  problem. 


4.2  Scaling  the  Model 

To  set  up  the  analysis,  we  must  prescribe  a  scaling  for  the  parameters  of  the  random 
medium,  in  this  case  the  permittivity,  permeabilit  y,  and  conductivity.  The  general  form 
is 

<(*)  ~  '"U)  = 

p(x)  =  pP(.r)  -  /<(. r.x/i')  (4.7) 

<t(x)  -  <t ( .r  )  rr{x,x/r) 

In  this  representation,  u  >  0  is  a  small  (<  .  •  1 )  dimensionless  parameter  which  charac¬ 

terizes  the  relative  scale  on  which  microscopic  variations  in  the  medium  take  place.  In 
previous  reports  we  have  shown  how  this  parameter  may  be  identified  in  terms  of  the 
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mean  free  path  between  elementary  scatterers  in  a  foliage  covered  area.4 

The  presence  of  two  scales  in  (4.7)  suggests  that  we  are  permitting  variations  of 
the  parameters  across  the  structure  of  two  different  spatial  types.  The  “fast  scale” 
y  =  x/u  captures  the  effects  of  microscopic  variations;  that  is,  as  we  pass  from  one 
elementary  scatterer  to  the  next.  Variations  on  the  “slow  scale”  x  capture  the  effects 
of  macroscopic  variations;  that  is,  variations  as  we  move  from  one  type  of  vegetation 
to  another,  or  from  one  region  of  a  given  type  to  another.  Note  that  as  y  changes  by 
one  unit,  the  macroscopic  scale  x  changes  by  v  <<  1  units;  and  as  r  changes  by  one 
unit,  y  changes  by  \/u  >>  1  units. 

Using  this  scaling,  there  are  two  basic  cases  which  can  be  treated  by  versions  of  the 
method  of  multiple  scales: 

1.  The  functions  e(sr,  y),  y.(x,  y),  a(x,  y)  are  periodic  in  their  second  argument. 

2.  The  functions  e(x,  y),  y.(x,  y),  a(r,  y)  are  random  functions  of  their  second  argu¬ 
ment  which  are  stationary  and  satisfy  a  mixing  condition. 

The  first  case  was  treated  in  (Bensoussan,  Lions,  and  Papanicolaou  1978)  and  in 
(Sanchez-Palencia  1980),  (see  also  the  work  of  Tartar  and  Murat  (Tartar  1986))  among 
many  other  related  problems.  The  second  case  provides  a  basis  for  modeling  certain 
types  of  foliage  covered  terrain. 

4That  is,  v  may  be  defined  in  terms  of  the  mean  free  path  between  scattering  centers  normalized  by 
a  characteristic  length  of  the  interaction  process,  e.g.,  the  length  of  a  scattering  region  or  the  principal 
axis  of  the  first  Fresnel  zone,  etc. 
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The  mixing  which  must  be  imposed  in  the  second  case  is  that  the  statistical  corre¬ 
lations  between  random  events  in  the  medium  must  decrease  rapidly  over  the  volume 
of  interest.  That  is,  the  correlation  between  values  of  t(x,y\)  and  t(x,y7)  decreases 
as  |yi  —  y2|  °o.  In  effect,  this  means  that  the  microstructure  of  the  foliage  at  one 

point  yj  cannot  depend  too  strongly  on  that  at  the  (distant)  point  y2.  This  is  clearly 
reasonable  in  many  situations.  It  may  not,  however,  be  me  case  in  wind  blown  grass. 

Notice  that  the  scattering  medium  is  still  permitted  to  have  (statistical  or  deter¬ 
ministic)  variations  which  take  place  in  the  “macroscopic"  x  spatial  scale. 


4.3  The  Multiple  Scale  Hypothesis 

The  objective  of  the  multiple  scale  asymptotic  analysis  is  to  compute  approximations  to 
the  fields  Ev(t,  x),  x)  in  the  limit  as  1/  — •  0,  that  is,  in  the  limit  as  the  microscopic 
variations  become  increasingly  dense. 

A  formal  procedure  for  computing  approximations  is  to  assume  an  expansion  of  the 
fields  in  the  form 

Ev{x,xlv)  =  £’0( x )  4-  i' £ i ( .r ,  y )  4  v2E7{x,y)  4-  •  •  ■  (4.8) 


with  y  —  x/v  and  a  similar  expansion  for  Hl '.  Introducing  the  change  of  coordinates 

(4.9) 


d 

dx 


d  ^  1  d 

dx  v  By 


in  Maxwell’s  equations,  substituting  the  expansions  for  Ev  and  Hv ,  one  finds  a  sequence 
of  equations  for  the  terms  Ekix,y)yk  ~  0.  1 ,2, .  .  .  (and  //*.).  These  may  be  solved  using 
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a  device  called  a  “corrector”  which  is  derived  from  a  kind  of  “separation  of  variables” 
argument. 


If  the  dissipation  is  zero  ( <r  =  0),  the  result  is  a  representation  of  the  solution  for 
the  approximate  field  in  the  form 

q<ix)  0  E0(x)  0  Vx  E0(t,x ) 

(4.10) 

0  q„(z)  H0{x )  J  —  Vx  0  H0(t,x) 

Eo(0,.r)  =  £o(*),jyo(0,r)  =  H0(x) 

V  •  [g«E0]  =  0,  V  •  [q^Ho]  =  0 


The  (matrix- valued)  functions  q((x),qfi(x),  which  depend  on  the  macroscopic  spatial 
scale,  are  defined  in  terms  of  c(x,  y), p(i,  y)  by  the  “corrected  averages” 

j  ,*,J  =  1,2,3  (4.11) 

where  is  either  or  and  the  function  \a(y)  is  the  corrector  associated  with 
e  or  p.  Note  that  the  effective  parameters  in  the  approximation  qt,q^  are  not  just 
the  averages  of  the  rapid  variations  of  the  parameters  e(.r,  •),  p(x,  •)  over  the  medium. 
They  include  the  “correctors”  which  retain  the  microscopic  interaction  effects  (multiple 
scattering)  in  the  final  approximation. 

The  correctors  satisfy  a  system  of  the  general  form  (see  the  following  sections  for 
details) 

V  •  [a(c_,  +  V \“(y))]  =  0,  A4[\,]  =  0  (4.12) 

where  ej  is  the  j,h  natural  basis  vector  in  S?3  In  these  equations  A4  is  the  operation  of 
averaging  over  a  typical  ceil  in  the  domain  for  the  periodic  case;  and  it  is  expectation 


q'J{x)  =  M  a, 


+ 


dx-(y) 

dy, 
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with  respect  to  the  stationary  distribution  describing  the  medium  in  the  random  case. 
When  the  conductivity  ir  /  0,  the  homogenized  model  is  more  complex,  including 
frequency  dependent  effects  -  see  below. 


4.4  Implementation  of  the  Approximation 

Thus,  to  implement  the  approximation,  we  must  solve  the  equations  (4.12)  for  the 
correctors  and  then  compute  the  averages  in  (4.11).  Computation  of  the  correctors  is 
possible  in  the  periodic  case  (see,  e.g.,  (Begis,  Duvaut,  and  Hassim  1981)  (Bougat  and 
Derieux  1978)).  However,  in  the  random  case  of  interest  here,  the  computation  is  very 
difficult.  The  equations  (4.12)  in  the  random  case  are  defined  path  by  path;  hence,  one 
would  have  to  have  a  complete  statistical  characterization  of  the  medium  to  be  able  to 
solve  them.  This  is  impossible  in  practice,  and  one  must  evaluate  the  correctors  by  a 
procedure  like  Monte  Carlo  simulations.  Hence,  these  equations  should  be  regarded  as 
the  basis  for  derivation  of  further  approximations  or  alternative  representations  for  the 
medium  which  require  only  physically  attainable  statistical  information,  e.g.,  second 
order  statistics  in  the  microstructure  of  the  medium.  Such  approximations  have  been 
derived  in  other  problems  (conductivity,  porosity,  etc.)  in  the  past.  We  shall  consider 
these  methods  and  their  adaptation  to  the  EM  scattering/absorption  problem  in  the 
next  section. 

Before  giving  the  details  of  the  derivation  of  the  approximation,  we  shall  make  a 
few  remarks  on  its  interpretation. 
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•  First,  note  that  the  general  form  of  the  approximation  is  valid  (with  a  very 
different  mathematical  interpretation)  in  both  the  random  and  periodic  cases. 


•  Second,  the  role  of  the  correctors  in  providing  the  “correct”  form  of  the  approxi¬ 
mation  cannot  be  dismissed.  Approximations  which  simply  average  the  variations 
of  the  permittivity  and  permeability  over  a  domain  are  incorrect  -  making  an  error 
of  “order  one”  in  the  small  parameter  v.  That  is,  the  approximation  (4.11)(4.12) 
is  an  approximation  of  “order  u"  (at  least  in  the  periodic  case);  omitting  the 
correctors  invalidates  this  estimate. 


4.5  Derivation  of  the  Effective  Parameters  Using  Homoge¬ 
nization 

There  are  two  basic  cases  to  consider: 

1.  The  periodic  case  when  the  inhomogenieties  in  the  medium  are  spatially  periodic; 
and 

2.  The  case  when  the  inhomogenieties  in  the  medium  are  randomly  distributed. 

The  expressions  for  the  effective  parameters  are  easier  to  grasp  in  the  periodic  case; 
and  they  are  more  readily  computable.  For  this  reason  we  shall  treat  this  case  first  as 
background  to  the  random  case. 
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4.5.1  Periodic  Media 


This  problem  was  treated  in  (Bensoussan  Lions  and  Papanicolaou  1978),  (Sanchez- 
Palencia  1980),  and  (Tartar  1979).  Our  treatment  is  based  on  (Sanchez-Palencia  1980). 


To  begin,  we  rewrite  Maxwell’s  equations 


dDv 

dt 


VxF-T  +  F 


dB u 
dt 


-V  x  Ev  +  G 


£w(z,0)=0,  Hv(x,  0)  =  0 


(4.13) 


where  F,G  are  (localized)  source  terms  and5 

£>r  =  Bf  =  ^»,-( ~ )Hj ,  Jf  =  <rii(l)Eur  (4.14) 

We  assume  that  €,p.,<r  are  smooth  Y  periodic  functions  of  x/u  which  are  symmetric 
and  positive  definite.6  Adopting  the  multiscale  hypothesis 


Ev(x,t)  =  E°(x,t)  +  vE\x,y,t)  +  ■  •  • 


(4.15) 


and  similarly  for  DV,B'J  and  Ju  with  each  term  in  the  asymptotic  expansion  periodic 
in  y.  We  shall  use  the  change  of  coordinates 


d  _  l_d 

dx  dx  ^  v  dy 


(4.16) 


5  We  use  the  summation  convention  throughout  this  section  -  repeated  indices  in  an  expression  are 
summed. 

GTliat  is,  =  e;t,  >  7l{p,7  >  0. 
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with  y  =  t/V,  and  similarly 


(V  X  •)  (Vr  x  •)  +  -  (Vy  x  •) 

v 

div  — *  divx  -f  -divv. 

v 

Rewriting  Maxwell’s  equations  in  terms  of  these  coordinates,  inserting  the  expan¬ 
sions  (4.15)  into  them,  and  equating  coefficients  of  like  powers  of  u,  we  arrive  at  a 
sequence  of  equations  for  E° ,  H° ,  Ex ,  Hl , The  0( i'-1)  and  CHI'0)  terms  are 

Vy  x  H°  =  0,Vv  x  E°  =  0  (4.17) 

8D° 

=  YX  x  H° +  Vyx  H'  -  J°  +  F  (4.18) 

at 

a  do 

w=-VrXli*-V,xB'+G 

The  equations  (4.17)  imply  that  Ho(*,y,t),  Eo(x ,  y,  t)  are  gradients  in  y  for  fixed  values 
of  the  other  arguments. 


The  terms  Vy  x  El  and  Vv  x  H 1  are  derivatives  of  periodic  (in  y)  functions;  hence, 
they  have  zero  mean  with  respect  to  the  averaging  operation 


M(Vy  x  El )  =  /,  Vy  x  Eldy  =  0 


(4.19) 


\y\  Jy 

(and  similarly  for  Vy  x  //’)  where  Y  is  a  typical  “cell”  of  the  periodic  structure  and 
|F|  is  the  volume  of  the  cell.  Therefore,  averaging  the  system  (4.18)  over  the  cell  Y , 
we  arrive  at  the  homogenized  system 

dD°  _ 

'dt 
8B° 


8t 


Vx  X  H°  -  J°  +  F 

=  -v,  x  e°  +  a 


(4.20) 
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where  the  overbar  denotes  the  averaging  operation.  To  complete  the  model,  we  must 
derive  the  homogenized  constituent  equations  corresponding  to  (4.14).  These  laws  take 
two  different  forms,  depending  on  the  presence  or  absence  of  dissipation  (<r  =  0  or  not). 


Consider  the  divergence  of  (4.13) 

div  +  7^  =  div  F 

div  =  aw  a 

Using  the  ansiaz  (4.15)  and  the  coordinates  (4.16),  we  have,  at  order  , 


i  r.( 


div 


div.  i  ®Ti  -  o 


and 

V  di 

which,  with  the  zero  initial  conditions,  implies 

divy  £?°  =  divy[//(x,  y)H°(r ,  y)]  =  0 

From  (4.17)  we  have  that  E°  and  H°  are  gradients;  hence,  we  can  write 


(4.21) 


(4.22) 


(4.23) 


(4.24) 


E°  -  E°  =  V<*>,  H°  -  H°  =  V t/< 


(4.25) 


with  the  means  of  <f>  and  V’  zero.  Using  these  in  (4.22)(4.23),  we  have 

dil' 


dy 


:  | Fijiy) 


jr~  + 

Pyj  } 


=  n 


d 

dy, 


+  ^j(y) 


dVJ  + 


=  0 


(4-26) 

(4.27) 
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These  are  the  key  equations  in  application  of  the  homogenization  method  to  Maxwell’s 
equations.  If  we  regard  (4.26)  as  an  equation  for  the  potential  V’,  then,  noting  that  H° 
does  not  depend  on  y,  we  can  solve  (4.26)  by  assuming  the  “separation  of  variables” 

0(*>! /,*)  =  \k{x,y,t)Hl{x,t)  (4.28) 

and  regarding  (4.26)  as  an  equation  for  \(y)  as  a  function  of  y,  treating  x,t  as  fixed 
parameters.  We  call  Xk(xi  l/>  Oi  ^  =  U  2, . . . ,  d  [d  is  the  dimension  of  the  space)  the  cor¬ 
rectors  associated  with  pij.  Using  (4.26)  and  (4.28),  the  equation  for  the  kth  corrector 
is 

d_ 

% 

This  problem  plays  a  key  role  in  the  homogenization  method.  We  call  it  the  cell  problem. 
The  assumption  that  p,j(x,y)  is  periodic  in  y,  symmetric,  and  positive  definite  (as  a 
matrix  for  each  x,y),  guarantees  that  (4.29)  has  a  unique  solution7  which  has  zero 
mean 


The  solvability  condition  for  (4,29)  is  that  the  right  hand  side  have  zero  mean.  This 
holds  since  p,j  is  periodic  and  smooth. 

The  homogenized  form  of  the  field  equations  may  be  obtained  using  the  correctors. 
Writing  out  the  0(v°)  term  in  the  expansion  of  (4.21),  we  have 

div^S1  +  divxB°  =  divxG  (4.30) 

7 Among  the  set  of  functions  periodic  in  y  with  first  derivatives  square  integrable  over  a  cell,  i.e. , 
periodic  functions  in  //|1<I,.(S?'1).  See  (Sanches-Palencia  1980),  pp.  51-54,  or  (Bensoussan  Lions  Papani¬ 
colaou  1978). 


/‘.j 


dxi 


dmk 

%  ’ 


k  —  1, 2, . . . ,  d. 


(4.29) 
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Using  M.  to  denote  the  averaging  operation,  the  solvability  condition  for  this  equation 
is 

^{div.fl0  -  div,G'}  =  0 

which  we  rewrite  as 

calM  [divse{^(x,y)[V!/tA  +  H°(x))  -  G}\  =  0  (4.31) 

Using  the  representation  (4.28)  for  the  potential,  (4.31)  holds,  if 


M  [MV^  +  DJ/o-C?]  =0. 


This  provides  the  definition  of  the  homogenized  magnetic  permeability,  fih,  as 


/4  =  M 


y )  + 


(4.32) 


The  effective  permittivity  and  conductivity  will  depend  on  the  frequency  of  the 
incident  radiation;  and  the  analysis  must  reflect  this.  Recall  (4.22),  rewritten  as 


div  | 

[r* 

+  crE0|  =  0 

(4.33) 

=  div  | 

—  6  +  a  [ V„<£  +  Euj  j 

=  ^divv  teV^]  +  divi/ 

(4.34) 

-t-divv 

f  9r 

+  divv  [(tE0] 

Thus,  the  local  “cell”  equation  is  an  evolution  (in  /)  equation. 
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To  treat  (4.34)  we  use  a  technique  introduced  in  (Sanchez-Palencia  and  Sanchez- 
Hubert  1978).  Let  Vy  be  the  space  of  suitably  smooth  functions  0(y)  which  are  periodic 
with  zero  mean,  and  let  Vy  be  endowed  with  the  inner  product 

f  d<j>  80  , 

<4>,6  > vY  —  /  €ij—  —  ay 

Jy  dyi  dyj 

Using  this,  (4.34)  becomes 

dr  T  d<f>  —  1  90 

di  Jy  [eiid^+eij{Eo)j\ 

r  d<f>  80  ,  r  .  -  .  dO 

+  +  Sr^Eo)^r° 

Using  the  inner  product  notation,  this  may  be  written  compactly  as 

d 


di 


<  <t>  +  Cj(Eo)j,0  >vv  +  <  A<j>,0  >vY  +  <  c2j(Eo)j,  0  >vY—  0 


where 


A  %  <Ji 


d 


'rU‘dy.iy 

do 


tJdyj 

<  Clj,0  >Vy=  jy 

<  c),0  >vy  =  jY  (Tij—dy 
The  evolution  equation  (4.36)  can  be  solved  for  the  potential  <f> 

m=  -c)(Eo)j(i)  +  J*  e~A^~T^c*(Eo)j(T )dr 


C*  =  4c! 


,  s2 


(4.35) 


(4.36) 


(4.37) 


Using  this  expression  for  <f>,  we  can  compute  the  average  flux  vector  D0  and  average 
induced  current  J0  and  define  the  effective  permittivity  and  conductivity.  Using  Eq  = 
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E0  4-  V<f>,  we  write  the  average  of  Do  =  ( E0  as 

(Do),  =  M{et}(Eo  +  V<*>)}  (4.38) 

=  M{ellt(E0)lt} 

~(-c)(E0):(t)  + J^e-^c^EoU^dr^  j 

This  may  be  written  compactly  as 

(Do),  =  a\j(E0)j{t)  +  b'j(t  -  T)(E0)j(t  -  T)dr  (4.39) 

where 

dc)(y) 

a*j  =  M  tijiy)  -  (4.40) 

Kj  =  M  («‘4r4)  (»)  (4.41) 

In  a  similar  fashion,  we  can  compute  the  induced  current 

(Jo),  =  a^EoW)  +  fb°3(t  -  t)(Eo)M  -  r)dr  (4.42) 

J  0 

where  a a  and  ba  are  defined  by  (4.41)(4.42)  with  c  replaced  by  <r. 

Using  Laplace  transforms,  we  can  more  readily  express  the  effective  scattering  pa¬ 
rameters  as  a  function  of  frequency.  Indeed,  if  we  let  Dl/( u>)  be  the  Laplace  transform 
of  Du  with  respect  to  f,  and  similarly  for  the  other  variables,  then  Maxwell’s  equations 
take  the  form 

wD"  +  r  =  V  v  +  F  (4.43) 

u iBv  =  -V  x  Ev  +  G 
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The  claim  is  that  Maxwell’s  equations  approach  the  homogenized  laws 

cuDh  4-  Jk  =  (cue  4-  <r)hEh  (4.44) 

Bh  =  iihHh 

where  the  superscript  (-)h  denotes  the  “homogenized”  laws.  Using  the  evolution  equa¬ 
tions  (4.39)(4.42)  we  can  identify 

cua-j  4-  cub‘j(u>)  4-  afj  4-  ^j(ou)  =  +  <7t;)\  ftcu  >  0  (4.45) 

This  completes  the  derivation  of  the  (first  order)  homogenization  theory  for  Maxwell’s 
equations. 

Remark:  Using  the  arguments  in  (Bensoussan  Lions  and  Papnicolaou  1978)  (Sanchez- 
Palencia  and  Sanchez-Hubert  1978)  and  (Sanchez-Palencia  1980),  it  can  be  shown  that 
the  homogenized  Maxwell’s  equations  (4.39)(4.42)  (or  the  frequency  domain  version) 
have  a  unique  solution.  Moreover,  it  can  be  shown  that  the  homogenized  Maxwell’s 
equations  are  the  limits  as  v  — ♦  0  of  the  original  scaled  Maxwell’s  equations  (in  an 
appropriate  weak*  topology). 

Rather  than  give  these  arguments,  we  shall  turn  our  attention  to  the  case  when  the 
coefficients  in  Maxwell’s  equations  are  random  processes  in  the  spatial  variables. 

4.5.2  Random  Media 

As  we  have  noted  earlier,  the  formulas  for  the  effective  parameters  in  the  homogenized 
form  of  Maxwell’s  equations  carry  over  in  form,  at  least,  to  the  random  case.  The 
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averaging  operation  A4  is  interpreted  as  expectation  with  respect  to  an  appropriate 
probability  measure.  While  the  procedure  of  simply  expanding  the  field  quantities  in 
asymptotic  series  in  the  small  parameter  is  only  formal  in  the  random  case,  it  can 
still  be  used  with  appropriate  cautions  to  determine  equations  for  the  homogenized 
field  quantities.  The  arguments  for  convergence  of  the  scaled  field  quantities  to  the 
homogenized  quantities  are,  however,  substantially  different  in  the  random  case. 

The  central  difficulty  in  treating  random  coefficients  arises  in  the  meaning  and 
evaluation  of  the  “correctors”  x',X'“,  X° •  They  are  random  processes  which  must  be 
evaluated  (numerically)  by  a  Monte  Carlo  type  procedure.  Before  numerical  methods 
are  designed  it  is  important  to  address  the  fundamental  issue  of  existence  and  unique¬ 
ness  of  solutions  to  the  corrector  equations.  Since  each  of  the  three  correctors  vVxM)  \ ° 
satisfies  the  same  type  of  elliptic  equation  with  random  coefficients,  we  shall  treat  the 
generic  problem:  Find  \(y;w)  such  that 


x(0;u;)=0,  E\{y)  =  0,  Vy 

dx/dyj  is  a  square  integrable  stationary  process  Vy  and 

d  (  ,dX\  dgj 

in  the  sense  of  distributions.  We  assume  that 


(4.46) 


ao|£l2  <  £  °ij(y,w )&•£,■  <  — 1£|2  G  9T,a0  >  0 
ij  30 


(4.47) 


where  atJ(y,u>)  and  gj{y,u> )  are  square  integrable  stationary  processes,  i,y  =  l,...n. 
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We  shall  consider  a  set  up  for  stationary  processes  as  presented  in  (Papanicolaou  and 
Varadhan  1979).  Let  (H,.4,  P)  be  a  probability  space  characterizing  the  uncertainties 
in  the  underlying  medium,  and  define  H  =  L2(ll,  A,  P).  We  assume  that 

there  exists  a  strongly  continuous  unitary  group  Ty  on  H,y  €  9?"  (4.48) 

Tv  is  ergodic,  which  means  if  /  £  H  satisfies  (4.49) 

Tyf  =  /,  Vy,  then  /  is  a  constant. 

if  /  >  0,  then  Tyf  >  0  and  Tv  1  =  1  (4.50) 

The  group  Ty,  which  is  used  to  define  a  stationary  processs  on  the  probability  space, 
has  a  spectral  resolution  defined  by 

Ty=  f  e'XvU(d\) 

J&n 

where  U(d\)  is  a  projection  valued  measure.  We  next  define 

A7M  =  (Tyf)  (u>)|v=0  (4.51) 

which  are  closed  densely  defined  linear  operators  with  domains  T>(Di)  in  Ti.  Let  Ji1  = 
C\U  ^ Di)  wWch  is  dense  in  "H.  We  equip  T~0  with  the  Hilbert  scalar  product 

((/,$))*'  =  Efg  +  Y,  EDJDj9  (4.52) 

i 

We  identify  H  with  its  dual  and  call  H~l  the  dual  of  Hl.  We  have  the  inclusions 

Hx  C  H  C 

each  space  being  dense  in  the  next  one  with  continuous  injection.  The  family  Ty  is  also 
a  strongly  continuous  unitary  group  on  Hl . 
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Remark  The  periodic  case. 


Let  Q  be  the  unit  n  dimensional  torus,  A  the  ^-algebra  of  Lebesgue  measurable  sets 
and  P  Lebesgue  measure  on  U.  Then  Ti  is  the  space  of  measurable  periodic  functions 
(period  1  in  each  component)  such  that 

/  ( f(uj))2<Lo  <  oo. 

Jn 

We  define 

Tyf{“)  =  f(*  +  y) 

hence 


Dj 


d 

d<”i 


/• 


An  important  fact  in  the  periodic  case  which  does  not  carry  over  to  the  random  rase, 
is  that  there  is  no  analogue  of  Poincare’s  inequality  (Bensoussan,  Lions,  and  Papani¬ 
colaou  1979). 


Consider  now  random  variables  not  necessarily  in  Ti.  We  assume  that 


Ty  is  a  linear  group  on  the  set  of  complex  random  variables  such  that 
Vt)]  . .  .fjk  complex  random  variables,  x['  Borel  bounded  function  on  Ck ,  then 

E^Tyiji,. . .  ,Tyf)k)  -  E'Mvi - - fjk)  (4.53) 

y,u>  — »  Tyf)  is  measurable, 

Tyij  >  0  if  ij  >  0. 
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A  ilationary  process  is  a  stochastic  process  T)(y;u} )  which  can  be  represented  in  the 
form 

Tj(y;u>)  =  TyTjiu).  (4.54) 

The  space  of  square  integrable  stationary  processes  can  be  identified  with  7i .  Moreover, 

=  DxTyf]{uj)  =  TvDif]  (4.55) 

if  rj  £  T~0  . 

Note  that  the  continuity  assumption  on  Ty  implies  that  the  square  integrable  stationary 
processes  are  necessarily  continuous  function  of  y  with  values  in  T~C.  Hence,  if  fj  £  H, 
n{y\u)  €  C,0(9?n;7i),  space  of  uniformly  continuous  functions  on  5?"  with  values  in  H. 
If  fj  £  then  rj(y\ u>)  £ 

The  Cell  Problem  We  consider  here  stationary  processes  atJ(j/;u>)  such  that 

a°|£|2  <  <  ~ -|£|2  V^e»",ao>0  (4.56) 

“  On 

Let  gAVi^)  =  Tyijj  square  integrable  stationary  processes,  j  =  l,...n.  We  shall  solve 
the  problem:  Find 


X(y;u>)  £  <->(»";?*), X(0;u>)  =  0,£*(y)  =  0  Vy 
d\/dyj  is  a  square  integrable  stationary  process  Vj 


(4.57) 


in  the  sense  of  distributions  with  values  in  7 ~i  (or  as  continuous  functions  with  values 
in  H~ 1 ). 


Papanicolaou  and  Varadhan  (1979)  have  shown  the  existence  and  uniqueness  of  the 
solution  of  (4.57).  We  shall  follow  their  proof  with  minor  changes.  Note  that  x(y;u>) 
itself  is  not  a  stationary  process.  This  is  a  big  difference  with  respect  to  the  periodic 
case  and  relates  to  the  Remark.  Note  also  that  x(y;u>)  €  C2(9?"; H-1).8 


Let  Xj  £  H  such  that 


we  can  assert  that 


(y;u>)  =  Ty\j 


ExjDk<f>  =  E\kDj4>,  Vtfre'H'. 


Indeed  we  have  to  show  that 


as  an  equality  in  H  1 .  But 


hence 


Dk\j  =  DjXk 


TyDkXj  =  Xj(y;w) 
vyk 


d2X 

dyjdyk 


(y;w) 


TyDkXj  -  TyDjXk 


which  implies  (4.59). 


(4.58) 


(4.59) 


We  have  also 


®By  virtue  of  (4.50)  it  is  sufficient  to  have  (4.5fi)  for  u  =  0. 
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Ex  =  0. 


(4.60) 


This  follows  from 

\(y\u))  =  52  f0  TfvXjyjM 

hence 

B(x/»)  =  E  =  o 

i 

by  the  assumption.  Therefore  (4.60)  follows. 

The  main  result  in  this  section  is  the  following: 

Theorem.  There  exists  one  and  only  one  solution  of  the  corrector  equation  (4.57). 

This  means  that  the  problem  of  computing  the  effective  parameter  for  the  homog¬ 
enization  problem  by  averaging  the  coefficient  aij(x,u>)  with  respect  to  the  underlying 
probability  measure  is  well-posed.  See  (Bensoussan  and  Blankenship  1988)  for  an 
application  of  this  result  in  homogenization  of  nonlinear,  stochastic  control  problems. 

Proof. 

Uniqueness. 

Assume  that  dgj/dyj  =  0.  Define 

(4.61) 

Note  that  (4.57)  can  be  written  as 

- D,Ty(alvXj)  ~  0  in  H~l 
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hence,  as  is  easily  seen 


EauxjDi^  =  0. 


(4.62) 


Because  of  (4.59)  we  have 


A  kV(d\)xj 

|iA  -  0\2 


\jU(d\)xk 

HA-/3P 


which  implies 


-  &  =  /„  £ 


(-i\i  -  l3)(i\j  -  f3) 


U(d\)xk(u>)  =  Xk- 


Therefore,  (4.62)  reads 


EavXjXi  +  tlEdyXj^  =  °- 


However, 


E\l3^\2  =  / 


X 


(-(A;  -  f3){-i\k  -  0)) 
|tA-/?|« 


EU(d\)xjXk 


-  £Ef/{0}.u\fc  as  /5  — >  0 
and  by  ergodicity  and  property  (4.60),  we  have 


E\0<t>a\2  -  0  as  0  -  0. 


Therefore,  (4.63)  implies 


Ed,j\jX,  =  0 


(4.63) 
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hence  x»  =  0,  which  implies  also  \  =  0. 


Existence: 


Let  (3  >  0,  we  solve  the  problem 


d  i  dx*  \  ,  a  0  d9> 
ay  dyo  oyj 

x‘‘ir^)  =  TvxB,x>  eW‘ 


This  problem  is  equivalent  to 

Ea,yDjx0Dt4>  +  l3Ex&4>  =  -EgjDj4>,  V<£  e  H' 


We  easily  deduce  the  estimates 


E|£>a"I!  <  C 


l3E{xBf  <  C. 


Let  us  extract  a  subsequence  such  that 

Dj\$  —i >  Xj  in  H  weakly. 
EDjX0Dk<j>  =  EDk\0  D)4> 

ED:x 0  =  0 

we  deduce  (4.59)  and  (4.60).  Going  to  the  limit  in  (4.65)  we  have 


EdijXjDj 4>  =  -Eg}Dj4>  V<£  £  Hl . 


(4.64) 


(4.65) 


(4.66) 
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Define  then 

X(y;t*>)  =  fx,Se'Xy  ~  -  i^j)V(dX)\j[uf)  (4.67) 

then 

=  Ty\k 

and  x(0;u>)  =  0,  Ex(y)  =  0.  Then  (4.66)  can  be  written  as 

—  D,a,jXj  =  D,g,  equality  in  H~l 

which  is  indeed  (4.57). 

In  (Bensoussan  and  Blankenship  1988)  this  result  was  used  to  study  the  homge- 
nization  of  nonlinear,  stochastic  control  problems.  In  (Begis  and  Blankenship  1987) 
numerical  methods  for  the  evaluation  of  homogenized  representations  to  certain  prob¬ 
lems  in  the  mechanics  of  period  composite  materials  were  discussed  (see  also  (Begis, 
Duvaut,  and  Hassim  1981)).  The  problem  of  designing  numerical  methods  for  the 
treatment  of  random  media  remains  open. 
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5  A  Macroscopic  Scattering  Model 


In  this  section  we  describe  a  computational  model  for  the  scattering  and  absorption  of 
EM  radiation  by  a  volume  of  foliage  covered  terrain.  In  effect  we  combine  the  detailed 
models  of  the  previous  sections  with  a  model  for  the  foliage  as  a  “very  rough  interface” 
between  the  air  and  the  underlying  earth.  The  effective  parameter  models  developed 
in  the  earlier  sections  are  used  to  characterize  the  scattering  and  absorption  of  the  EM 
radiation  within  the  microstructure  of  the  foliage.  The  result  is  an  effective  parameter 
representation  for  the  scattering  volume.  A  rough  surface  (interface)  scattering  model 
is  then  used  to  characterize  the  interaction  of  the  radiation  with  the  “layer”  of  foliage. 
We  use  the  method  of  smoothing  perturbation  (Keller  1964)  to  treat  this  “macroscopic” 
aspect  of  the  interaction. 

Typical  regions  of  foliage  covered  terrain,  viewed  as  a  scattering  surface,  are  very 
rough  in  the  sense  that  they  do  not  satisfy  the  validity  conditions  of  the  classical 
Kirchhoff  integral  method  and  the  Rayleigh-Rice  perturbation  methods.  In  particular, 
scattering  regimes  in  which  the  rms  value  h  of  the  surface  height  variations  is  substan¬ 
tially  larger  than  their  correlation  length  £,  while  A,  the  wavelength  of  the  incident 
radiation  is  comparable  to  h  (“the  low  frequency  case”)  or  to  f  (“the  high  frequency 
case”)  are  incompatible  with  the  classical  methods.  We  shall  call  these  surfaces  “very 
rough.”  The  local  angles  of  incidence  between  very  rough  surfaces  and  the  incom¬ 
ing  radiation  are  not  always  small,  even  in  operational  geometries  which  suggest  near 
“grazing”  incidence.  Hence,  it  is  necessary  to  have  methods  suitable  for  all  angles  of 
incidence. 
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Our  solution  technique,  which  is  an  approximate  method,  involves  several  steps. 
First,  the  foliage  covered  terrain  is  regarded  as  a  rough  surface,  which  is  in  turn  regarded 
as  constituting  a  “layer”  of  inhomogeneous  material  separating  two  infinite  slabs  of 
(essentially)  homogeneous  material  (e.g.,  air  and  earth).  Using  this  paradigm,  the 
problem  is  reformulated  as  scattering  in  a  random  medium  in  which  the  inhomogeneity 
is  characterized  by  a  small,  dimensionless  number  u  =  \/h ,9 

The  method  of  smoothing  perturbation  is  applied  to  the  vector  wave  equation  for 
the  electric  field  vector  in  the  entire  space  in  a  manner  similar  to  the  work  of  Keller 
and  Karal  (1964).  This  leads  to  a  deterministic  wave  equation  for  the  average  field 
which  involves  the  second  order  statistics  -  the  correlation  function  -  of  the  surface 
fluctuations.  Solution  of  this  equation  may  be  readily  carried  out  using  reflection 
and  transmission  coefficients  for  the  random  layer,  and  these  parameters  constitute  an 
effective,  simplified  representation  of  the  original  scattering  problem  which  requires  a 
minimum  of  empirical  data  for  implementation.  In  numerical  studies  we  have  used  the 
technique  to  treat  surfaces  with  uniform  and  Gaussian  distributed  height  variations.10 
The  reflection  coefficients  differ  only  slightly  from  the  Fresnel  coefficients  at  grazing 
incidence;  however,  at  larger  incident  angles  substantial  differences  (20-50%)  in  both 
amplitude  and  phase  variations  are  obtained. 

9This  parameter  is  in  addition  to  the  small  parameters  characteriiing  mean  free  path  length  in  the 
random  medium  models  used  in  the  previous  sections. 

1(,See  (Blankenship  1986)  and  (Barakati  1983).  A  sample  of  these  results  is  included  below. 
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5.1  The  Method  of  Smoothing  Perturbation 

The  mathematical  technique  we  use  is  a  standard  method  in  perturbation  theory  which 
has  a  simple,  general  formulation  (Keller  1964).  Let  u o  be  a  wave  in  a  homogeneous 
medium  characterized  by  a  linear  operator  Lo ,  so  LoUo  =  0,  with  suitable  boundary 
conditions.  Let  (Q ,F,P)  be  a  probability  space,  and  let  u>  6  fl  designate  a  different 
medium  characterized  by  the  operator  Lo  -  i'Li(aj)  -  ^2L2(^)  4-  0(is3).  Here  v  is  a 
measure  of  the  inhomogeneity  of  the  medium,  and  L\  and  Li  are  perturbation  operators 
representing  the  effects  of  inliomogeneity.  A  (random)  wave  in  this  medium  satisfies 

[L0  -  vL\(u>)  -  i^2L2(u>)  +  O^Ju'V)  =  0.  (5.1) 

Let  <  u  >  be  the  expected  value  of  u( u>)  as  a  random  variable  on  (Cl,F,P), 

<  u  >=  f  u(y)P{du>). 

J  n 

Our  objective  is  to  find  an  equation  for  <  uv  >  accurate  to  0( v3).  When  <  L\  >=  0, 
we  can  write 

[Lo  -  u\<  L,L0- ‘L,  >  +  <  Li  >)]  <  uv  >=  0{u 3)  (5.2) 

When  the  0{vz)  terms  are  omitted,  this  is  an  explicit,  deterministic  equation  for  <  uv  > 
which  is  accurate  to  0(v3)  (see  the  remarks  on  convergence  below).  This  equation  may 
be  treated  by  a  number  of  methods. 

In  many  applications  uv(x, u>)  is  a  vector  function  of  the  (vector)  position  i.  Then 
Lo,  Lj,  and  L2  are  matrix  operators,  each  component  of  which  may  be  a  differential 
operator.  In  these  cases  Lq1  is  an  integral  operator,  the  kernel  of  which  is  the  (dyadic) 
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Green’s  function  (?o(x,y),  a  matrix  function  defined  by 

L0G0(x,y)  =  I6{t  -  y)  (5.3) 

where  I  is  the  identity  matrix  and  6  is  the  Dirac  delta  function.  Then,  in  terms  of  the 
Green’s  function  we  have 

L0{x)  <  uv{x)  >  - V 2  <  Tj(a-)  J  Go(x,y)Li(y)  <  uu(y)  >  dy  >=  0{v3)  (5.4) 

It  is  clear  from  this  how  the  second  order  statistics  -  cross  correlations  and  autocorre¬ 
lations  -  enter  the  express  r  <  v1'  >. 

Notice  that  the  smoothing  perturbation  leads  to  an  integral  equation  for  the  ap¬ 
proximation  of  <  u  >.  This  equation  is  more  complex  than  the  differential  equations 
which  arise  in  other  approximation  methods.  This  complexity  leads  to  greater  accuracy 
and  to  an  increased  range  of  validity  in  the  approximation.  The  papers  (Chow  1974) 
(Blankenship  1979)  contain  a  comparative  analysis  of  several  approximations,  including 
the  smoothing  perturbation. 

5.2  Random  Surface  Model 

Consider  the  physics  of  a  plane,  time-harmonic  EM  wave  of  arbitrary  polarization 
incident  on  a  very  rough  interface  separating  air  and  a  lossy  dielectric  (foliage).  The 
rough  interface  is  represented  by 

2  =  /(r^)  =  (5-5) 

v 

r  -  ( .r ,  y ) 
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where  g(r,u)  is  a  spatially  homogeneous  random  field,  and  u>  £  fl,  ( fl,F,P )  a  proba¬ 
bility  space,  is  a  realization  of  that  random  field.  We  assume  that  j<7(r,u>)|  <  hj 2  for 
some  h  >  0  and  Vr.  The  dimensionless  parameter  v  >  0  is  the  correlation  length  of 
g(r,tj})  normalized  by  h  .  We  are  interested  in  the  surface  scattering  problem  as  u  — +  0. 
We  shall  assume  that  the  random  surface  has  zero  mean 

<S(r, •)>=().  (5.6) 

This  puts  us  in  the  framework  corresponding  to  equations  (5. 2, 5.4).  Spatial  homogene¬ 
ity  means  that 

<  g(r<-)g(s,-)  >=  R(r  -  s).  (5.7) 

And  if 

(  =  [f  R(r)dr/R(  0)j1/2 

JR* 

is  the  correlation  length  of  the  surface  variations  (recall  h.  =  H(0)),  then  u  =  (/h.  The 
example 

R(r)  =  (5.8) 

_  h2e^  |r‘IA') 

with  r  =  r/h  illustrates  the  definition. 

The  surface  cross  section  shown  in  Figure  5.2  illustrates  the  scattering  geometry 
and  the  parameters.  For  a  fixed  h,  as  v  — *  0,  the  surface  fluctuates  more  rapidly  about 
the  plane  z  —  0,  and  the  incident  radiation  encounters  an  increasingly  rugged  interface. 
As  the  figure  suggests,  the  surface  variations  are  confined  to  a  “layer”  of  thickness  h. 
We  shall  use  this  artificially  defined  layer  as  the  basis  for  our  analysis. 
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Figure  5.1:  Hough  surface  scattering. 

5.3  Vector  Wave  Equation  and  the  Smoothing  Approximation 


Starting  with  Maxwell’s  equation  for  time-harmonic  fields,  we  can  obtain  the  following 
vector  wave  equation  for  the  electric  field  vector  E 

V  x  V  x  Ev  -  k2Nl(r,u>)El/  =  0  (5.9) 


valid  in  the  entire  space  r  =  (x,y,z).  Here  k  =  2n/\  is  the  propagation  constant  in 
free  space,  and  7Vj(r,  u;)  is  the  random  dielectric  constant  associated  with  the  layered 
structure  defined  in  Figure  5  2.  We  have  assumed  that  the  magnetic  permeability  p  is 
uniform  (p .  =  1)  throughout  space.  Defining  the  function 


1  z>gv(r,  u>) 

< 

0  elsewhere 


(5.10) 
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then 

Nl{r^)  =  /T(r, u,)  +  r?2[  1  -  fT(r, u,)\  (5.11) 

where  tj2  =  er  j60\cr,  a  the  conductivity,  is  the  complex  dielectric  constant  of  the 
material  in  region  2  of  Figure  5.2.11 

The  single  (vector)  equation  (5.9)  together  with  a  radiation  boundary  condition  for 
Eu(r)  as  |r|  — »  oo  is  a  reformulation  of  the  EM  interface  problem  which  includes  a  wave 
equation  for  z  >  gl/(r,u>),  a  wave  equation  for  z  <  gt/(r,<jj),  and  the  usual  continuity 
relation  for  the  fields  in  regions  1  and  2  across  the  interface.  Thus,  (5.9)-(5.11)  is  a 
reformulation  of  the  Maxwell’s  equation  with  a  random  boundary  condition  as  a  partial 
differential  equation  with  a  random  coefficient.12  This  formulation  enables  us  to  apply 
the  smoothing  perturbation  io  (5.9).  Let 

Hq{z)  =<  Hv(r,-)  > 

N20(z)  =  H0(z)  +  q2[l  -  H0(z)\  (5.12) 

*/l/(r,u>)  =  (JEF'(r,u;)  —  H0(z)\/v 

(recall  the  assumption  of  spatial  homogeneity).  Then  (5.9)  may  be  rewritten  as 

VxVxT-  k2N^(z)Eu  -  u2k2{\  -  i}2)Hv{r,uj)E''  =  0.  (5.13) 

nThe  complex  dielectric  constant  would  be  computed  or  estimated  using  the  methods  discussed  in 
the  previous  sections. 

12Note  the  similarity  to  the  models  used  in  the  section  on  homogenisation.  In  effect,  the  models  here 
are  “large  scale"1  representations  of  the  scattering  phenomena;  whereas  those  used  in  the  homogenisation 
technique  are  local,  “small  scale"1  representations.  It  would  be  possible  to  combine  the  rough  surface 
paradigm  with  the  model  used  to  represent  microscopic  events;  however,  the  resulting  analysis  would  be 
rather  complex. 
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Identifying 

L0  =  Vx  S/~k2N2(z) 

Lx  =  ~k2{  1  —  if  )Hl’(  r,  u; ) 

L2  =  0 

and  noting  <  Lx  >— -  0,  we  can  apply  the  smoothing  perturbation  to  (5.13).  That  is, 
<  Ev  >  satisfies,  formally  to  0( ix3 ), 


Lq  <  Ev  >  -v2  <  I.I-'L,  ><  Ev  >  =  0 


(5.15) 


Let  G0(r|f)  be  the  dyadic  Green’s  function  of  the  operator  L0  with  an  outgoing  wave 
condition  (Tai  1971);  i.e., 

L0Go(r\f)  =  I6(r\f)  (5.16) 


We  can  then  write  (5.15)  as 

L0  <  Ev  >  (r)  -  u2k\  1  -  if)2  j  j  I  G’0( 

X  -  X  y-y 


r  r , 


where 


■R 


R 


V  V 


x  -  x  y-y 


\  dx  dy 

i)  <  Ev  >  (r)  —  ~dz  =  0 


v  v 


:<  //>,•)//>,)  > 


v  v 

Setting  x  —  (x  -  x)/u.y  -  (y  -  t/)/^,5  --  z,  we  can  rewrite  (5.17)  as 


L0  <  E1'  >  (r)  -  ,/V(l  r?2'2 


if3 


G0(r/.r,  i'y,z,z) 


(5. 17) 


(5.18) 


(5.19) 


■  R(x .  y,  z.  z)  <  El  >  (x  --  ux .y  -  vy.  z  )dxdydz  —  0 
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Now  assume  that  R  is  integrate  over  (x,y)  £  ??2.  (See  Appendix  I  for  comments  on 
this  assumption.)  Then  for  <  E1'  >  bounded,  the  integral  in  (5.19)  wili  be  bounded 
since  iZ(x,y,  z,£)  —  0  for  |r|,  |£|  <  h/2. 

Adopting  the  ansatz 

<  Ev  >  =  <  E  >0  W  <  E  >2  <  E  >4  +  . . .  (5.20) 

and  substituting  in  (5.19),  we  find 

Lo<E>q=0  (5.21) 

U  <  E  >2=fc4(l  -t)2)2  I"12  Go(0,0,z,z)  /  f  R(x,y,  z,  z)dxdy  (5.22) 

J-h/2  J  JW 

•  <  E  >0  (x,y,z,z)dz 

Thus,  <  E  >0  approximates  the  average  electric  field,  formally,  to  0(v2).  An  0(n3) 
approximation  may  be  obtained  by  substituting  <  E  >0  into  (5-22)  and  solving  for 
<  E  >2-  We  shall  not  consider  this  extension  here. 

5.4  Reflection  Coefficients 

To  illustrate  the  computations  involved  and  the  qualitative  effects  of  the  surface  rough¬ 
ness  on  reflectivity,  we  shall  compute  the  reflection  coefficients  for  a  surface  whose 
deviations  from  the  mean  plane  are  uniformly  distributed  in  [  —  h/2,  h/2}.  This  means 
that  H0{z)  is  linear  in  \ -h/2,  h/2 and,  in  effect,  it  presumes  minimum  information 
about  the  detailed  structure  of  the  rough  surface.  We  shall  treat  horizontally  and  ver¬ 
tically  polarized  incident  waves  separately.  Arbitrary  polarized  incident  waves  may  he 
treated  bv  superposition. 


Let  the  waves  be  incident  in  the  x  -  z  plane  and  let  9  be  the  angle  of  incidence  of  the 
(plane)  wave.  (See  Figure  5.4.)  Let  5,  and  v,  denote  the  unit  vectors  in  the  direction 
of  the  incident  electric  field  vector  for  Horizontal  and  vertical  polarization,  respectively. 
Similarly,  hr  and  vr  are  the  corresponding  unit  vectors  for  the  reflected  field.  Let  k, 
and  kr  denote  the  unit  vectors  along  the  directions  of  propagation  of  the  incident  and 
reflected  waves,  respectively.  Using  i,y  and  z  to  denote  the  unit  coordinate  vectors, 
we  can  define  the  polarization  by 

h,  -  hT  ~  x  x  i  =  -y 

Vi  =  k  x  ki,  >\  =  K  x  5r  (5.23) 

In  the  case  being  considered 

k,  -  sin  Ox  -  cos  Oz,  kr  -  sin  Ox  cos  Oz  (5.24) 

and  this  in  (5.23)  means  that 

{>,  —  cos  Ox  -  sin  Oz,  x\  —  cos  Ox  -  sin  Oz  (5.25) 


’.4.1  Horizontal  Polarization 


In  this  case  the  average  electric  field  vector  in  c  >5/2  can  be  written  as  a  combination 
of  the  incident  and  reflected  waves  as 


E 


■  n  - 


jtr\T  5in  0  (  r 


h  i  2 } 


r  °'9]h, 


^  (pH  H  h , 


.  flr’j  sill  ^  ;  h  /2|«"OS0i 


>  5/2 


7(j 


(5.26) 


where  pnH  and  puv  are  the  reflection  coefficients  associated  with  horizontally  polarized 
incident  radiation.  The  coefficient  puv  will  be  zero  if  there  is  no  depolarization  due  to 
scattering.  (This  is  the  case,  at  least  to  order  zero  in  v.) 

Having  assumed  uniformly  distributed  surface  deviations,  the  vector  wave  equa¬ 
tion  (5.21)  for  the  zeroth  order  average  scattered  field  is 

VxVx<£>oTV(;)<£>o=0  (5.27) 

where 

1  z  >  h/2 

f(z)  =  l^-z  +  1^2-  jz|  <  h/2  (5.28) 

r?2  z  <  -h/2 

From  this  we  see  that  (5.26)  is,  in  fact,  a  solution  of  (5.27,5.28). 

The  field  in  z  <  —h/2  can  be  written  as 

<  E  >0=  (7HJ/^  +  ^HVt\)eifcl'co,'-u+fc/2)v^-,ina#,1  r  <  -h/2  (5.29) 

where  we  have  required  <  E  >0  to  be  outgoing  and  finite  as  z  — *  —  oo.  The  unit  vectors 
ht  and  vt  denote  the  transmission  directions  of  the  horizontally  and  vertically  polarized 
waves,  respectively.  They  are  defined  by 

h  z  =  -y 

k,  =  A- sin  Or  -  A  -'Jr;7  -  sin2  0:  (5.30) 

kt  x  ht  \/q2  -  sin2  <9  sin  0 
Ifc/i  >]  V 

The  constants  *)hh  and  ~i  hv  are  transmission  coefficients. 
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The  electric  field  in  |;|  <  hj 2  can  be  written  as  a  linear  combination  of  two  vector 
wave  functions,  M  and  N  (Tai  1971) 


M  =  V  x  (V’i) 


(5.31) 


N  = 


kf(z) 


V  x  V  x  (</>i) 


which  satisfy  (5.27,5.28)  if  we  take  the  scalar  functions  V’  and  <p  as  solutions  of 


V2tH  fc2/(*)V'  -  0  (5.32) 


2  ,  1  df(  z  dip  2 

v  <t>  -  7—  —T—  /( =  0 

f(z)  dz  dz 


Assuming  a  separation  of  variables 


</•  -  ^(^(y^U) 


<t>  =  X2(x)Y2(y)Z2(z) 

and  separation  constants  k2a2  and  k2b2 ,  we  obtain  from  (5.32) 


(  —  +*V  )*,(*)  =0,i  =  1,2 
(~+k7b2)Y,(y)  =  0,r  ----  1,2 

[  +  k2{f{z)  -  a2  -  b2)\Zx{z)  =  0 

azl 

1  df(z)  d  _  2  2  1.2 , 


ldz2  /(;)  d;  dr 
The  form  of  the  incident  wave  suggests 


~  +k2(f(z)-a2  -b2))Z2(z)  =  0 


(5.34) 


(5.36) 

(5.37) 


sin  6-0 


(5.38) 
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so  that 


.Y,(t)  :=  i  --  1 , 2 

Y,{y)  -  1,  i=I,2 

The  general  solution  of  (5.36)  may  he  written  using  special  functions  as 

Zi(z)  =  W1/3Mff;;>(u.)  +  BH[%(w))  (5.40) 

where 

2 

u<  =  w(s)  -  —{cz  +  d)3/2  (5.41) 

3  c 

c  =  k2  ,  rf=-t2(l  f  r?2 )  -  A.-2  sin2  0 

h  2 

and  =  1,2,  are  the  Hankel  functions  of  the  i-th  kind  of  order  p  (Abramowitz 

and  Stegun  1971 ). 

The  second  order  equation  (5.37)  with  f(z)  linear  cannot  be  reduced  to  a  standard 
form.  We  can,  however,  construct  two  linearly  independent,  series  solutions  Z2i(z i) 
and  Z22{zy)  in  terms  of  Cj  --  f{z),  Jzj  <  hj 2;  see  Appendix  II  for  details.  Note  that 
the  values  of  Z2\  and  Z22  are  needed  at  only  two  points,  which  makes  the  numerical 
computation  particularly  attractive. 


Writing  a  general  solution  of  (5.37)  as 


Z2  ~('Z2\(z] ) D Z22{ z j  ) 


[5.42) 


we  can  write  a  general  expression  for  the  average  electric  field  in  the  layer  |r|  <  hj 2 

c,  r.  .  ^  (  1  -  'f  K';2l(‘l  )  4  Dz22(Zi) 

<  £  ()-  |j  sin  t9—7 - - - — — 

h  /( ^ ) 
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-  \jkan0wl,3{AH[ya[w)  +  (5.43) 

+  \kfnn*9(CZn{zi)  +  DZ22{zx))z\,\z\  <  h/2 

where  A,B,C and  D  are  complex  constants  to  be  determined  from  the  boundary 
conditions  at  — -  ±h/2. 

The  boundary  conditions  satisfied  by  <  E  >0  are  derived  from  the  requirements 
that  the  tangential  components  of  the  electric  and  magnetic  fields  must  be  continuous 
across  the  (artificial)  “boundaries”  of  the  “layer”  at  z  —  h/2.  These  conditions  reduce 
to  the  requirements: 


(i)  zx  <  E  > o  is  continuous  at  z  —  h/2  and  z  —  —h/2 

(ii)  2  X  Vx  <  E  >0  is  continuous  at  r  =  h/2  and  2  =  -h/2. 


Each  of  these  conditions  yields  four  equations,  so  we  have  eight  linear  equations 
relating  the  eight  unknowns  Phh,  Phv,  A,  B,  (7,  D.  ~ihh,  and  7 hv-  These  equations  de¬ 
couple  into  two  sets  of  equations,  and  they  may  be  written  as 


-1 

«12 

0- 13 

0 

PHH 

1 

0 

®22 

023 

-1 

A 

0 

]  k  cos  8 

032 

033 

0 

B 

-jkc  os  8 

0 

042 

043  jk 

yji]2  -  sin2  8 

7  HH 

0 

(5.44) 


where 

a,2  =  «’i./3^!/3(«’+)>  Q13  =  U>+/3HJ/’(U>+) 

022  ~  U>1/3//J^(ll-_  ).  023  =  U’1/3H|/3(u»_) 
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a32  =  (d  +  ch/2)'l2wl^H%3{w,  ),  o33  =  (d  +  ch/2)'/2w\/3H{]}/3(w+  ) 
ai2  =  (d  -  ch/2)1/2w\/3H("/3(w _  ),  o43  -  (d  -  ch/2)l/iwlJ*H%9{w-  ) 


—  1  b\2  b\ 

0  b22  b2 

— ik  cos  6  b32  b3 

0  642  bi 


(5.45) 


it i7  -  sin2  9  ~,nv 


where 


612  —  Z2  i(l),  5]  3  —  Z22O) 

b22  =  -Z2i(rj2),  b2  3  =  ~Z22(t]2) 

V  V 

h  (1  m  a  (1  ~  ***)£  m 

O32  —  7  22l(  1  ),  f>33  ~  T  “22 (  1  ) 


1  U  -  »?’).  ,  2.  ,  (1  -  »?2)-  ,  2. 

642  —  —  -21  ( 27  ),  043  —  ^  ^22(2/  ) 


and  u?±  =  w(z  =  ±h/2).  In  writing  (5.44,5.45)  we  have  replaced  A  and  B  by  jkAsinS 
and  jlcB  cos  0,  respectively.  We  have  also  absorbed  a  factor,  ksh\9,  into  C  and  D  and 
used  the  fact  (Abramowitz  and  Stegun  1971)  that 


W'3H\'hw))  =  ,(«•),  .  =  1,2 


(5.46) 


and  dw/dz  =  ( cz  -f  d)2 .  The  roots  of  complex  quantities  are  interpreted  as  principal 


values.  For  convenience  we  write  (5.44,5.45)  in  compact  notation  as 


PhQhh  =  * 


P\  < Ihv  -  0 


(5.47) 


(5.48) 
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with  the  definitions  of  Ph,Pv >  etc.  obvious  from  (5.44,5.45). 

Numerical  results  show  the  matrix  Py  in  (5.45)  to  be  nonsingular.  So  we  obtain 

Phv  =  C- '  =  D  =  -)hv  —  0  (5.49) 

for  horizontal  polarization.  This  means  that,  to  zero  order  in  the  small  parameter  u, 
there  is  no  change  in  polarization  in  scattering  or  transmission.  The  reflection  coefficient 
pUH  is  computed  from  (5.44). 


5.4.2  Vertical  Polarization 

The  field  in  z  >  h/2  may  be  written 

<  E  >0=  eJM*  «n*-(i-fi/2)  co»«]  (5.50) 

+  (pHvhr  +  p V V  f'r  ) €",fe|*  # +< * ~ h/ 2 )  C°* ^ ,  Z  >  h/2 

The  expression  for  the  field  in  |z|  <  h/2  is  unchanged  from  (5.43);  and  in  z  <  —  h/2  we 
have 

<  E  >0=  (7 vaht  +  z  <  -h/2.  (5.51) 

Applying  the  boundary  conditions  at  z  —  ±j,  we  obtain 

PaqvH  =  0  (5.52) 

Pyqvv  ~  s  (5.53) 
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where 


qvH 


PVH 

pvv 

A 

C 

B 

i  qw  — 

D 

7vh 

ivv 

(5.54) 


with  Pff,Pv,  and  5  as  in  (5.44,5.45,5.47,5.48).  Once  again  numerical  computations 
show  that  there  is  no  depolarization  to  zero  ordei  in  the  small  parameter  v  =  l fh\  that 
is,  pyu  =  A  —  B  =  7vh  =  0.  The  reflection  coefficients  may  be  computed  from  the 
4x4  system  (5.53). 


5.5  Typical  Numerical  Results 

Figures  5. 5-5. 5  show  the  magnitude  and  phase  of  phh  and  pyy  computed  from  (5.47) 
and  (5.53),  respectively,  for  a  very  rough  interface  with  er  =  80,  cr  =  4.3mho/m,  A  = 
2.5cm  and  h  =  A/2,A,2A.  For  comparison  we  have  included  the  magnitudes  and  phase 
of  the  complex  Fresnel  reflection  coefficients  for  horizontal  and  vertical  polarization 
(Northam  1981 ). 


cos  8  —  \j i)2  ~  sin2  8 

PFH  - - ■■■■-■=  ■■■== 

cos  8  +  y  i]2  —  sin2  8 


(5.55) 


Pfy 


i]2  cos  8  -  J ?/2  -  sin2  8 


T]2  cos  8  +  \Jt)2  -  sin2  8 
The  plots  show  the  coefficients  against  the  grazing  angle,  90°  —  8,  rather  than  the  angle 

of  incidence  8. 


lb)  hi  =  2.5  c  m. 
(c)  h=  5.0  cm. 


_J _ i _ 1 _ 

30°  60° 

grazing  angle 


>•  I:  ('Iihm'  <>|  i>hh  jirnzinj>  nuujp  ami 


grazing  angle 

Figure  5.5:  Magnitude  of  pvv  versus  grazing  angle  and  h. 

Evidently,  the  differences  between  the  reflection  coefficients  of  a  very  rough  surface 
and  the  Fresnel  coefficients  of  a  flat  surface  become  pronounced  as  the  grazing  angle 
increases.  At  normal  incidence  the  magnitude  of  the  reflection  coefficient  \phh\  of  the 
very  rough  surface  is  about  20%  larger  than  that,  \pFH\,  of  a  flat  surface.  The  phase 
difference,  Figure  5.5,  is  even  more  pronounced.  And  the  differences  are  exaggerated 
even  further  in  the  case  of  vertical  polarization.  We  have  shown  the  dependence  of  phh 
and  pvv  on  h  for  h  =  A/2,  A,  2 A.  We  find  that  as  h  decreases  the  magnitude  of  the 
reflection  coefficient  decreases  and  the  phase  increases.  The  Fresnel  coefficients  define 
the  limits  as  h  — ►  0. 
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I'  if*ur«-  Yfi:  1'liasr  >>l  I 


\  e.r;r/inn  an 


S  X 


At  very  low  grazing  angles  (less  than  5°)  the  magnitude  and  phases  of  the  rough 
surface  reflection  coefficients  differ  little  (<  3%)  from  the  corresponding  quantities  for  a 
flat  surface.  However,  the  very  rough  surface  scattering  model  is  useful  in  representing 
regions  of  very  rough  geometry  on  a  large-scale  surface  which  otherwise  satisfies  the 
requirements  of  the  Kirchhoff  method.  In  this  case  the  local  angle  of  incidence  at  the 
very  rough  surface  regions  will  assume  all  values  in  the  range  0°  — 90°.  Over  this  range  of 
angles  the  reflection  coefficients  derived  from  the  very  rough  surface  scattering  model 
depart  significantly  from  the  Fresnel  reflection  coefficient.  Ignoring  these  differences 
can  lead  to  significant  signal  processing  errors. 


5.6  Reflection  Coefficients  for  a  Gaussian  Surface 

For  non-uniformly  distributed  surface  profiles  f(z)  is  not  linear;  but  we  can  construct 
a  piecewise  linear  approximation  to  the  actual  f(z)  and  then  extend  the  methods  of  the 
preceding  paragraphs  to  the  resulting  “multilayer”  scattering  problem.  For  example, 
suppose  the  surface  profile  g(r,u>)  is  Gaussian  with  zero  mean  and  standard  deviation 
s.  We  set 

U  =  3s  (5.56) 

which  is  based  on  the  fact  that  the  probability  of  a  Gaussian  random  variable  taking 
values  outside  (—3s,  3s)  is  less  than  1%. 

Next  we  construct  a  piecewise  linear  approximation  to  f{z)  the  Gaussian  probability 
distribution  function.  Since  f(z)  is  an  odd  function  of  z,  we  need  an  odd  number  of 
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linear  segments  to  approximate  the  profile.  The  approximation  used  is 


7772(2  +  h/2)  if  —  h/2  <  2  <  -3 

+  s)  4-  m2{h/ 2  -3)  if  —  s  <  z  <  3  (5.57) 

7712(2  —  3)  -f  27/7 1 3  +  m2(h/2  -  3)  if  s  <  2  <  h/2 

1  1  -  v/'2  /j» 

mi=v^’  m2=_ il 

In  this  case  we  have  approximated  the  layer  in  which  surface  variations  occur  (see 
Figure  5.2)  by  three  layers  each  with  a  linear  profile.  Writing  out  the  field  expressions 
in  terms  of  reflection  and  transmission  coefficients  and  free  parameters  ( A,B,C ,  etc.) 
and  applying  the  boundary  conditions  (i)  (ii)  at  the  (four)  artificial  boundaries  leads 
to  a  total  of  16  equations  for  the  coefficients  of  an  arbitrarily  polarized  wave.  For  a 
wave  with  either  horizontal  or  vertical  polarization  we  have  8  (non-trivial)  equations. 
In  general,  for  an  n -segment  piecewise  linear  approximation,  one  must  solve  two  sets 
of  (2n  +  2)  equations  to  compute  psHH  and  ply-  The  numerical  implementation  is 
straightforward. 

In  Figures  5.6  and  5.6  we  show  the  reflection  coefficients  and  for  horizontally  or 
vertically  polarized  radiation  incident  on  a  surface  with  a  Gaussian  profile  approxi¬ 
mated  by  three  linear  sections.  The  relevant  parameters  are  3  =  0.2  cm,/?  =  1.2cm,  A  = 
2.5cm,  cr  =  80,  and  <r  --  A.Zmho/m.  W'e  have  also  included  plots  of  the  Fresnel  coef¬ 
ficients  and  the  reflection  coefficients  for  a  linear  profile  distribution  with  h  ~  1.2cm. 
Note  that  the  reflection  coefficients  are  almost  identical  at  low  grazing  angles  where 
the  scattering  physics  are  insensitive  to  the  surface  profile.  The  coefficients  differ  sig¬ 
nificantly  at  larger  grazing  angles.  This  means  that  in  situations  like  remote  sensing  of 


the  ocean  surface  or  in  high  altitude  radar  engagements  or  in  cases  involving  composite 
surface  models  where  the  (local)  grazing  angles  could  take  values  well  away  from  0,  a 
careful  description  of  the  surface  roughness  distribution  is  important. 


5.7  Applicability  of  the  smoothing  perturbation  model  to  a 
Gaussian  surface 

Our  purpose  here  is  to  provide  a  formal  argument  supporting  use  of  the  smoothing 
perturbation  model  for  modeling  random  rough  surfaces.  We  verify  that  the  method 
applies  to  a  Gaussian  surface  satisfying  appropriate  conditions.  For  simplicity,  we 
consider  a  one-dimensional  rough  surface,  z  =  g(x,u>). 

Earlier  we  mentioned  that  the  smoothing  perturbation  model  applies  to  a  random, 
very  rough  surface  which  is  represented  by  a  zero-mean,  strict-sense  stationary  random 
process,  defined  on  a  probability  triple  (fi,  F,P),  provided 

f  R(x ,  2 ,  z)dx  <  oo  (5.58) 

J  —  oo 

where 

R(x,z,z)  =<  H(x,z,u)H(x  -  i,z,u>)  >  (5.59) 

H(x,z,u>)  =  H(x,z,u,)-H0(z)  (5.60) 

with  H{x,z,u>)  defined  to  be  one  for  z  >  g(x, u>)  and  zero  otherwise.  In  (5.60),  f/0(z) 
denotes  the  expected  value  of  H{x,z,u>). 
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PHASE  MAGNITUDE 


I'igniT  •  >.>v  M iiunit mlc  ami  phase  of  />\  \  for  Gaussian  surface. 


From  (5.59)  we  can  write  R(x,z,z)  as 


R(x,z,z)  =  <  H(x,z,u>)H(x,x,z,u>)  >  -H0(z)H0{z) 

-  P{g(x,“>)  <  z,g(x  -  *,u>)  <  S)}  -  H0(z)H0(z)  (5.61) 

=  F7{z,S,z)-H0(z)Ho{z) 

where  F7(',  ■;  x)  denotes  the  joint  distribution  of<7(r,u>)  and  g(x  —  i  ,w).  By  definition, 
Ho(-)  is  the  distribution  function  of  g(x,u). 

Suppose  the  autocorrelation  properties  of  the  random  process,  g(x,u;),  are  such 
that  the  random  variables  g(x,u> )  and  g(x  —  x,u;)  are  uncorrelated  for  x  >  d  and  for  all 
x  £  S?.  Then  we  have 


Fi(z,z,x)  =  H0{z)H0{z),  |*|  >  d.  (5.62) 

Using  (5.61)  and  (5.62),  we  obtain 

R(i,z,z)~  0,  |*|  >  d  (5.63) 

in  which  case  (5.58)  is  satisfied. 

It  is  reasonable  to  expect  that  (5.58)  will  also  be  satisfied  for  surfaces,  g{x, u>)  with 
an  autocorrelation  function 


Rg{*)  =<  g(x,u>)g(x  -  *,w)  >  (5.64) 

which  is  finite  for  all  f  and  decays  to  zero  as  |*|  —»  oo.  Since  R(x,z,z)  given  by  (5.61 ) 
depends  on  the  joint  distribution  function  of  the  rough  surface,  we  can  provide  only  a 
formal  argument  supporting  this  contention. 
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When  g{x, u))  is  Gaussian,  we  can  explicitly  write  the  joint  distribution  function, 
F2{-,-;x)  of  g(x,u>)  and  g(x  -  r,u>)  as 


F2(z,z,;x)  =  f  f  - 7 ===^dudv 

J -00  J-oc  2nJs*  —  R2(i) 


~?nv(  r)uf /*'  +  t 


(5.65) 


where  s  =  /?p(0)  denotes  the  variance  of  the  Gaussian  surface  roughness  g(x,uj).  Us¬ 


ing  (5.65),  we  can  write  (5.61)  as 


-r  c-^L.u1+erJ^ 

J~co  v  2ns  2  l  ,/2(s4 


.3(2  -  Rg(x)v/s2) 


f2 ~  Rl(i)) 


~Ho(z)Ho(z) 


where 


1  2  -t1  j 

7=e  1  dt 


(5.66) 


(5.67) 


Now  suppose,  the  autocorrelation  function  of  the  surface  roughness,  Rg(x),  satisfies 


Rg(x)  <  s2e  for  some  a  >  0. 


(5.68) 


Given  a  small  6  >  0,  we  can  then  find  L  such  that  for  lil  >  L 


Rgi*) 


<  6e-ali-Lj. 


(5.69) 


Inspecting  (5.66),  we  see  that  f^LR(x,z,z)  =  K(L )  will  be  finite.  Using  (5.69),  we 
find  that  for  |x|  >  L ,  the  argument  of  the  error  function  in  (5.66)  is 


s(z  -  RAx))v/s2  1  /  ,  ,  \ 

v  L-  =  ■—  (1  -  tl  3(2  _  iWe'"|r~£|)) 

12 (s4  -  Rl(x )  v/2.9  V  / 


v/2.9 


-q|x-£.| 


+  0(6z)e 


2 
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so  that  by  expanding,  the  error  function,  erf ,  in  a  Taylor  series  at  we  get 

r  \  26j>e-23/2’Vo|i~t| 


..A* -**(*) |’/»!) .. 

fii?  -  flj(i) 


7T  3 


(5.70) 


+  0(62)e 


2  v„-q|*-£| 


Using  (5.70)  we  can  write  (5.58)  as 

I*0  R(x,z,z)dx  =  K(L)  +  /  rfi(  f  - 
J-OO  J\r\>L  7-0 C  y/2w 


i  C~v  '2‘2  if  z 

dv  -  {  1  +  erf  — —  - 
its  2  \  \/2 a 


(5.71) 


r  _  /  2  ^2  - 

-J5To(«)J5ro(*)  -  r-e~o|*~L|  f  ve~"3/is2dv  4-  0(tf *)e-°,*-r',J 

ZtT  S  *  •'“Oo 


Using 


in  (5.71 ),  we  get 


r*  e‘'uJ/2,J  1  / 

ffo(i) = L  =  5  (' + cr/ 

f°°  R(x,z,z)  =  A'(I)  +  0(«)  <  c 
J-00 


y/2  s 


(5.72) 


(5.73) 


which  shows  that  the  smoothing  perturbation  model  applies  to  Gaussian  very  rough 
surfaces  with  autocorrelation  functions  satisfying  (5.68). 
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6  Conclusions 


The  purpose  of  this  work  was  to  develop  analytical  tools  for  the  modeling  and  analysis 
of  scattering  of  EM  radiation  in  regions  of  foliage  covered  terrain.  We  have  focussed  on 
developing  desciptions  of  the  multiple  scattering  which  occurs  within  a  region  contain¬ 
ing  a  dense  array  of  elementary  scatterers.  The  three  methods  we  have  examined  for 
this  problem  each  have  strong  points:  The  T-matrix  method  leads  naturally  to  approx¬ 
imations  for  the  scattered  field  based  on  an  explicit  representation  of  the  interaction  of 
elementary  scatterers  (of  several  different  classes).  The  computation  of  bounds  based 
on  the  geometry  of  the  underlying  material  provides  tests  for  the  evaluation  of  the 
approximate  representations  -  e.g.,  approximations  which  provide  estimates  near  the 
upper  or  lower  bound  on  the  effective  parameter  model  of  the  region  suggest  that  the 
model  is  valid.  The  homogenization  technique,  which  is  valid  for  low  frequency  radi¬ 
ation,  provides  not  only  a  systematic  method  for  the  construction  of  approximations 
to  the  scattering  physics;  it  also  provides  analytical  methods  for  the  evaluation  of  the 
precision  of  those  estimates.  That  is,  it  is  possible,  at  least  in  principle  to  estimate 
the  quality  of  the  approximation  based  on  rigorous  convergence  arguments.  Numerical 
evaluation  of  the  models  computed  by  homogenization  for  periodic  media  are  straight 
forward.  However,  the  case  of  random  media  is  quite  different  and  much  more  difficult. 

We  have  also  developed  a  computational  procedure  for  computing  the  EM  field 
scattered  by  a  “very  rough”  random  surface  (above  a  region  containing  a  dense  array 
of  elementary  scatterers).  In  this  scattering  regime,  the  rms  value  h  of  the  surface 
height  variations  is  substantially  larger  than  their  correlation  length  f,  while  A,  the 
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wavelength  of  the  incident  radiation  is  comparable  to  h.  The  method  of  smoothing 
perturbation  is  used  to  obtain  a  deterministic  vector  wave  equation  for  the  average 
electric  field  scattered  by  the  random  “very  rough”  surface.  The  coherent  scattering 
is  characterized  by  reflection  coefficients  which  can  be  computed  easily.  The  technique 
is  applied  to  treat  surfaces  with  uniform  and  Gaussian  distributed  height  variations, 
and  the  results  are  compared  to  the  Fresnel  coefficients.  Only  small  differences  are 
observed  at  grazing  incidence;  however,  at  larger  incident  angles  substantial  differences 
(20-50%)  in  both  amplitude  and  phase  variations  are  obtained.  Since  the  interaction  of 
EM  radiation  with  such  surfaces  involves  local  angles  of  incidence  within  the  full  range 
(0  —  90°),  errors  can  occur  in  approximations  which  ignore  these  differences. 
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